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Summary 


This  final  report  summarizes  the  research  carried  out  on  the  development  of  numerical 
approaches  for  simulating  the  plumes  of  pulsed  plasma  thrusters.  A  significant  amount  of 
progress  has  been  made  in  the  three  years  of  this  grant  and  this  is  summarized  in  the  9  journal 
and  conference  papers  that  are  included  as  appendices  to  this  report.  Our  modeling  has  made 
progress  in  all  aspects  of  simulating  these  complex  devices  including  Teflon  ablation,  plasma 
formation,  electro-magnetic  acceleration,  plume  expansion,  and  particulate  transport.  Several 
different  pulsed  plasma  thruster  devices  were  modeled  including  an  electro-thermal  device  (the 
PPT-4  developed  at  the  University  of  Illinois  by  Dr.  Rod  Burton),  a  Z-pinch  device  (the  AZ-PPT 
developed  at  Princeton  University  by  Dr.  Eddie  Choueiri),  and  various  micro-PPT's  (developed 
at  the  Air  Force  Research  Laboratory  by  Dr.  Greg  Spanjers  and  colleagues). 
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Electrical  Discharge  in  the  Teflon  Cavity  of  a  Coaxial 

Pulsed  Plasma  Thruster 

Michael  Keidar,  Member,  IEEE,  Iain  D.  Boyd,  and  Isak  I.  Beilis,  Member,  IEEE 


Abstract — In  this  work,  we  analyze  the  physical  processes  of 
a  pulsed  discharge  in  a  dielectric  (Teflon)  cavity.  This  type  of 
discharge  is  generated  in  a  coaxial  pulsed  plasma  thruster  (PPT) 
having  a  central  Teflon  cavity  to  produce  a  high-pressure  cloud 
of  ablation  products  during  the  discharge  pulse.  The  primary 
intended  role  of  this  model  is  to  provide  upstream  boundary 
conditions  for  particle  simulation  codes  used  to  study  the  exhaust 
plume.  The  main  features  of  the  electrical  discharge  in  the 
dielectric  cavity  include  Joule  heating  of  the  plasma,  heat  transfer 
to  the  dielectric,  decomposition  of  the  dielectric  followed  by 
partial  ionization,  and  acceleration  of  the  plasma  up  to  the  sound 
speed  at  the  cavity  exit.  We  consider  a  diffuse  type  of  discharge 
assuming  that  all  plasma  parameters  are  uniform  in  the  cavity. 
The  system  of  equations  is  based  on  the  plasma  energy  balance, 
thermal  conductivity,  dielectric  ablation,  and  mass  balance.  It  is 
found  that  most  of  the  energy  of  the  plasma  column  is  carried 
off  by  particle  convection  to  the  dielectric  and  by  radiation.  It  is 
found  that  during  the  pulse,  the  electron  density  peaks  at  about 
10^^  m“®  and  decreases  to  10^^  m”®  toward  the  end  of  the 
pulse,  whereas  the  electron  temperature  peaks  at  about  2.2  eV 
and  decays  to  1.5  eV.  Teflon  surface  temperature  peaks  at  about 
650  K.  Predicted  plasma  temperature  and  ablated  mass  are  found 
to  be  in  agreement  with  available  experimental  data. 

Index  Terms — Ablation  controlled  discharge,  near-wall  sheath, 
pulsed  plasma  thruster  (PPT). 


1.  Introduction 

PULSED  plasma  thrusters  (PPT’s)  have  been  investigated 
since  the  early  1960’s  and  were  among  the  earliest  electric 
propulsion  systems  to  be  flown  by  the  United  States.  First, 
plasma  thruster  designs  were  based  on  the  plasma  accelerators 
that  were  developed  for  high-energy  plasma  injection  into 
thermonuclear  reactors.  Different  configurations  of  plasma 
accelerators  were  proposed  with  electromagnetic  [1],  [2]  and 
electrothermal  [3]  dominant  acceleration  mechanisms.  The 
ratio  between  electromagnetic  and  electrothermal  mechanisms 
of  acceleration  depends  on  thruster  geometry  and  electrical 
discharge  parameters.  Both  parallel  rail  electrode  and  coaxial 
concepts  were  developed  [4]-[6]  to  produce  the  thrust.  The 
principal  advantage  of  the  PPT  is  the  simple  propellant  system 
design,  which  provides  high  reliability.  However,  the  PPT  has 
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poor  performance  characteristics.  For  instance,  a  flight-qual¬ 
ified  PPT  design  had  an  efficiency  of  about  8%  [7].  Several 
directions  for  improvement  of  PPT  performance  are  under 
consideration  [8]. 

A  great  interest  in  the  development  of  small  satellites  occurs, 
which  causes  the  PPT  to  be  reconsidered  as  an  attractive  propul¬ 
sion  option  [9],  [10].  PPT’s  are  expected  to  provide  exact  im¬ 
pulse  bits  to  be  used  for  accurate  attitude  control.  PPT's  can 
provide  high  specific  impulse  of  1000  s,  impulse  bit  of  about 
30  /iN-s  per  Joule,  and  can  operate  with  arbitrary  low  power. 

In  order  to  assure  successful  operation  of  a  PPT  on  the  space¬ 
craft,  a  complete  assessment  of  the  spacecraft  integration  effects 
is  needed.  The  solid-fed  PPT  plume  contains  various  ion  and 
neutral  species  because  of  propellant  decomposition  and  pos¬ 
sible  electrode  erosion.  The  main  integration  issue  is  the  depo¬ 
sition  of  highly  condensible  PPT  plumes  on  spacecraft  surfaces. 
Some  attempts  at  PPT  plume  simulation  using  hybrid  direct  sim¬ 
ulation  Monte-Carlo  particle-in-cell  (DSMC-PIC)  simulations 
were  performed  recently  [1 1],  [12].  In  [1 1],  however,  some  arti¬ 
ficial  starting  conditions  were  employed.  Accurate  plume  mod¬ 
eling  requires  the  formulation  of  boundary  and  initial  condi¬ 
tions,  which  depend  on  the  specific  PPT  design  and  pulse  pa¬ 
rameters.  Thus,  it  is  necessary  to  analyze  the  physical  processes 
involved  in  plasma  generation  and  acceleration. 

In  this  investigation,  we  will  concentrate  on  a  specific  type  of 
PPT,  recently  developed  at  the  University  of  Illinois,  so-called 
PPT-4  [9],  [10].  This  PPT  has  a  coaxial  configuration  in  which 
Teflon  is  ablated  from  a  cylindrical  cavity  sitting  in  front  of  the 
central  electrode.  The  device  has  a  pulse  length  of  about  10  /xs, 
and  the  overall  specific  impulse  was  measured  to  be  830  s.  The 
main  physical  processes  in  this  type  of  PPT  occur  in  the  Teflon 
cavity.  Rapid  heating  of  a  thin  dielectric  surface  layer  leads  to 
decomposition  of  the  material  of  the  wall.  As  a  result  of  heating, 
decomposition  and  partial  ionization  of  the  decomposition  prod¬ 
ucts,  the  total  number  of  particles  increases  in  the  cavity.  The 
PPT-4  also  has  a  ceramic  nozzle  in  which  the  plasma  is  pre¬ 
sumed  to  be  accelerated  by  both  electrothermal  and  electromag¬ 
netic  effects  [9].  In  principal,  the  ratio  between  electromagnetic 
and  electrothermal  acceleration  mechanisms  can  be  changed  by 
using  different  operational  parameters,  such  as  cavity  length, 
pulse  form,  and  duration. 

In  the  present  work,  a  model  of  the  physical  process  of  plasma 
generation  in  the  dielectric  cavity  will  be  developed.  The  present 
model  allows  calculation  of  the  electron  temperature,  electron 
and  neutral  densities,  dielectric  surface  temperature,  and  the 
sheath  potential  drops  near  the  anode  and  dielectric.  The  present 
model  helps  us  to  understand  plasma  generation  and  energy  bal¬ 
ance  in  electrothermal  PPT’s.  These  calculation  results  will  be 
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Fig.  1.  Schematic  of  the  problem  geometry  and  eneigy  balance. 


used  in  other  plasma  plume  model  as  time-dependent  boundary 
conditions. 

n.  The  Model 

The  model  presented  here  describes  the  physical  processes 
of  pulsed  electrical  discharge  in  a  Teflon  cavity  sitting  in  front 
of  the  central  electrode-anode,  as  shown  in  Fig.  1.  The  main 
features  of  electrical  discharge  in  the  dielectric  cavity  include 
heating  of  the  plasma,  heat  transfer  to  the  dielectric,  decompo¬ 
sition  of  the  dielectric  followed  by  partial  ionization,  and  elec¬ 
trothermal  acceleration  of  the  plasma  up  to  the  sound  speed  at 
the  cavity  exit.  As  a  result  of  Teflon  decomposition,  the  plasma 
density  increases.  Thus,  the  stored  energy  is  expended  in  plasma 
generation,  plasma,  electrodes,  and  insulator  heating,  and  ac¬ 
celeration.  In  the  PPT-4,  the  two  electrodes  are  connected  with 
a  30®  half  angle  nozzle,  as  shown  in  Fig.  1  (the  annular  elec¬ 
trode-cathode  is  not  shown).  Because  of  the  specific  geometry 
of  this  thruster,  the  current  flow  is  mostly  parallel  to  the  plasma 
flow.  In  the  initial  stage  of  the  discharge,  the  interelectrode  gap 
is  filled  with  a  certain  amount  of  plasma  by  a  spark  plug.  The 
voltage-  and  surface-induced  variations  in  the  igniter  plasma 
may  lead  to  shot-to-shot  variation  of  PPT  thrust.  However,  usu¬ 
ally  the  energy  and  mass  of  the  first  ignition  plasma  constitutes 
no  more  than  a  few  percent  of  the  total  mass  and  energy  in¬ 
volved  in  the  discharge.  Therefore,  one  may  expect  that  the  exact 
conditions  of  the  igniter  plasma  will  not  substantially  affect  the 
performance  of  the  PPT  [9].  Because  no  experimental  data  for 
the  initial  plasma  density  exist,  in  the  present  model,  the  initial 
plasma  density  will  be  used  as  a  free  parameter  of  the  problem 
(see  Section  IB). 

The  main  features  of  the  discharge  in  the  PPT  Teflon  cavity 
are  analogous  to  the  capillary  discharge  [13]  and  ablation  con¬ 
trolled  arcs  [14]-[16].  During  the  discharge  pulse,  the  plasma  is 
heated  by  Joule  heat  and  is  cooled  by  radiation,  energy  losses 
from  particle  convection  to  the  anode  and  dielectric,  and  ioniza¬ 
tion.  Mechanisms  of  energy  transfer  from  the  plasma  column  to 


the  wall  of  the  Teflon  cavity  are  also  possible:  heat  transfer  by 
particle  fluxes  and  radiation  heat  transfer.  As  a  result  of  this  in¬ 
teraction  of  the  plasma  column  with  the  walls  of  the  cavity,  the 
walls  are  heated  rapidly.  The  temperature  of  the  surface  layer 
increases  and  reaches  the  critical  temperature  of  phase  transi¬ 
tion.  Further  analyzes  should  include  treatment  of  the  kinetics 
of  chemical  reaction  of  material  decomposition  [17],  [18].  The 
principal  mechanism  of  thermal  degradation  of  the  dielectric  is 
breaking  the  bonds  in  the  backbone  of  the  chain  [17],  which 
consumes  the  main  portion  of  energy  transferred. 

During  the  discharge  pulse,  the  nonuniformity  in  the 
discharge  distribution  across  the  Teflon  surface  may  cause 
overheating  followed  by  phase  change  of  the  propellant, 
which  causes  high  pressure.  High  plasma  pressure  in  the  PPT 
channel  may  lead  to  Teflon  macroparticle  generation.  However, 
the  mechanism  of  such  macroparticle  generation  from  the 
propellant  is  not  understood  sufficiently  for  modeling.  Another 
possible  source  of  macroparticles  is  the  spot  attachment  at 
the  electrodes.  Spot  attachment  at  electrodes  is  a  typical 
phenomenon  in  the  early  stages  of  discharge  [9],  [19].  The 
macroparticle  generation  phenomenon  is  beyond  the  scope  of 
the  present  study. 

During  the  discharge  pulse,  parameters  change  rapidly,  so 
an  important  problem  concerns  the  possibility  of  establishing 
local  thermodynamic  equilibrium  (LTE).  The  typical  PPT-4 
pulse  discharge  duration  is  about  10  pts.  In  a  homogeneous 
transient  plasma  [20],  complete  LTE  may  be  obtained  in 
0.3  jis  for  a  helium  plasma  with  electron  density  of  10^^  m”^. 
An  estimation  of  the  characteristic  times  for  ionization  and 
recombination  has  shown  that  the  ionization  and  recombination 
time  scales  for  ground  states  of  C  and  F  are  less  than  the  typical 
time  for  discharge  parameter  changes  (few  microseconds)  [21]. 
Therefore,  LTE  establishing  may  be  considered  during  the 
discharge  pulse.  The  calculation  of  relaxation  time  for  higher 
excited  levels  shows,  however,  that  the  quasisteady  state  is  not 
achieved  for  the  second  stage  of  ionization  [20].  An  estimate 
of  the  relaxation  time  for  elastic  collisions  has  shown  that 
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in  a  plasma  with  density  of  and  an  electron 

temperature  of  1-3  eV,  equilibrium  of  electrons,  ions,  and 
neutrals  is  established  on  a  time  scale  of  microsecond  [3]. 
Thus,  the  present  model  considers  temperature  equilibrium 
Te  ~Ti  =  T,  where  T  is  the  plasma  temperature. 

In  the  present  model,  we  will  assume  that  a  diffuse  discharge 
covers  the  entire  surface  of  the  anode.  The  anode  current  is  car¬ 
ried  by  electrons  and  controlled  by  the  potential  drop  in  the 
anode  sheath.  The  potential  drop  in  the  sheath  depends  on  the 
plasma  and  current  density  and  may  have  spatial  variation.  Gen¬ 
erally,  the  sheath  is  slightly  negative,  to  repel  the  excess  of  the 
thermal  electron  current,  so  that  the  electron  current  to  the  anode 
is  equal  to  the  circuit  current 

Ua  =  -T  Hhnll)  (1) 

where  T  is  the  plasma  temperature  in  [eV],  Uh  is  the  random 
electron  current,  and  I  is  the  circuit  current.  The  random  elec¬ 
tron  current  may  be  calculated  as 

hh  =  |ene(8eT/7rme)^/^7rH^  (2) 

where 

Tie  is  the  electron  density  at  the  anode  sheath-plasma  in¬ 
terface; 

rUe  is  the  electron  mass; 

e  is  the  electron  charge; 

Ra  is  the  anode  radius. 

The  electrostatic  sheath  also  appears  near  the  dielectric  sur¬ 
face  where  quasineutrality  condition  breaks  down.  The  transient 
sheath  formation  near  the  dielectric  surface  is  determined  by 
the  dielectric  permittivity.  In  the  considered  range  of  electron 
density  m~^),  the  characteristic  charging  time  [22] 

is  less  than  10”^°  s,  which  is  much  smaller  than  the  charac¬ 
teristic  time  of  the  discharge  parameter  changes.  Therefore,  we 
will  use  a  quasisteady-state  sheath  model.  Furthermore,  we  will 
assume  that  the  plasma  density  is  so  high  that  the  sheath  thick¬ 
ness  (^Debye  length)  is  much  smaller  than  the  electron  Larmor 
radius  (in  the  self  magnetic  field).  This  assumption  will  be  justi¬ 
fied  later.  As  in  the  case  of  the  anode  sheath,  the  sheath  voltage 
drop  Ud  should  be  negative  to  repel  the  excess  of  the  thermal 
electron  current,  so  that  the  electron  current  is  equal  to  the  ion 
current 


3)  The  plasma  column  is  in  local  thermodynamic  equilib¬ 
rium  (LTE). 

As  a  result  of  the  foregoing  assumptions,  the  following  gov¬ 
erning  equation  can  be  formulated  to  describe  the  time-depen¬ 
dent  behavior  of  different  plasma  parameters  and  the  relation¬ 
ship  between  them. 

The  electron  energy  balance  in  the  quasineutral  plasma 
column  reads 

J  n,  ^  ^  -  AZfn,niT^I\l  +  Xg) 

2  at  a 

-L{2T+Ua)-^{2T+U^  +  T)  (5) 

L/  iCa 

where  the  first  term  on  the  right-hand  side  represents  the  Joule 
heating  of  the  plasma  column,  j  is  the  current  density  j  = 
{I / IT Rl),  a  is  the  plasma  conductivity,  and  jd  is  the  ion  cur¬ 
rent  density  in  the  sheath  near  the  dielectric  jd  =  Ii/{27tRaL). 
The  second  term  describes  the  radiative  emission  energy  losses 
[23],  [24],  A  is  the  constant  (1.6  •  10^  in  SI  units),  n,  is  the 
ion  density,  and  Xg  =  Eg /Te  is  the  Eg  is  the  energy  of  the  low 
excited  state.  The  third  term  represents  the  energy  losses  caused 
by  electron  convection  to  the  anode.  The  last  term  represents 
energy  losses  caused  by  convection  of  electrons  and  ions  to  the 
dielectric. 

In  general,  the  thermal  regime  of  the  Teflon  will  be  con¬ 
sidered  to  consist  of  two  stages:  the  first  stage  involves  ini¬ 
tial  heating,  i.e.,  the  Teflon  surface  heating  up  to  the  onset  of 
phase  transition,  and  the  second  stage  of  Teflon  decomposition. 
A  third  stage  also  exists  called  late  ablation  [25],  i.e.,  Teflon  de¬ 
composition  after  the  pulse,  which  is  also  under  consideration. 

The  temperature  can  be  calculated  from  the  heat  transfer 
equation 

dT/dt  =  ad^Tjdx^  (6) 

where  a  is  the  thermal  diffusivity,  a  =  XjCpp,  where  A  is 
the  thermal  conductivity,  Cpis  the  specific  heat  capacity,  and  p 
is  the  specific  weight.  Equation  (6)  is  subject  to  the  following 
boundary  conditions  [18] 

^XdTldx{x  =  0)  -  q{t)  ^AHT-  Cp{Z  -  To)T 
XdT/dx{x  =:  oo)  =  0 

T{t^0)=To  (7) 


Ud  =  -Tln{I,j,/Ii)  (3) 

where  U  is  the  ion  current  given  by  the  expression 

li  =  0.25eZjne(8eT/7rma)^^^27ri?a^-  (4) 

Here,  Zi  is  the  charge  number,  m*  is  the  ion  mass,  and  L  is  the 
cavity  length  (see  Fig.  1). 

Starting  from  the  above  considerations,  a  simplified  model  of 
the  pulsed  discharge  in  the  dielectric  cavity  is  proposed  using  the 
following  basic  assumptions  (some  of  these  assumptions  will  be 
discussed  in  Section  IV). 

1)  We  have  considered  a  diffuse  type  of  discharge  assuming 
that  all  plasma  parameters  are  uniform  within  the  plasma 
column. 

2)  The  plasma  is  quasi-neutral. 


where 

a:  =  0  corresponds  to  the  inner  dielectric  surface; 

AH  is  the  ablation  heat; 

r  is  the  ablated  flux; 

To  is  the  initial  room  temperature; 

q{t)  is  the  density  of  the  heat  flux; 

consisting  of  the  radiative  and  particle  convection  fluxes,  deter¬ 
mined  according  to  (5),  and  Tg  is  the  Teflon  surface  temperature. 

In  the  first  stage,  which  according  to  calculation  is  about 
1.5  ps,  the  energy  losses  connected  with  decomposition  are 
small  and,  thus,  may  be  neglected  in  the  energy  balance.  We  will 
use  an  approximate  solution  for  temperature  within  the  Teflon 
[17],  [18],  which  at  the  Teflon  surface  reads  [21] 


T,{t)^To  +  ^ 


12A 


W)CpP 


f\{t) 

Jo 


dt 


0.5 


(8) 
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After  the  surface  reaches  some  critical  temperature,  material  de¬ 
composition  begins  and  ablation  heat  becomes  significant  in  the 
energy  balance.  The  principal  mechanism  of  thermal  degrada¬ 
tion  of  the  Teflon  is  the  breaking  the  bonds  in  the  backbone  of 
the  chain  [17].  Heat  transfer  analysis  shows  that  the  temperature 
profile  for  the  ablated  Teflon  is  exponential  [17] 

T{x)  -  Ts  exp{-xTCp/X).  (9) 

We  will  use  this  approximate  solution  to  estimate  the  tempera¬ 
ture  gradient  at  the  Teflon  surface.  It  should  be  noted,  however, 
that  the  calculation  shows  little  influence  of  the  temperature  gra¬ 
dient  on  the  Teflon  surface  temperature. 

Using  (9)  with  boundary  conditions  (7),  it  is  possible  to  cal¬ 
culate  the  temperature  of  the  Teflon  surface  for  known  ablated 
flux  r,  which  in  turn  depends  on  the  equilibrium  pressure.  For 
known  surface  temperature,  one  can  calculate  the  equilibrium 
pressure  using  the  Teflon  formula  [9] 

P  =  P,  exp(-T,/T,)  (10) 

where  P  is  the  equilibrium  pressure  and  Pc  and  Tc  are  the  char¬ 
acteristic  pressure  and  temperature,  respectively.  In  the  quasis¬ 
teady  state,  the  equilibrium  pressure  should  be  equal  to  the  total 
pressure  in  the  discharge  column 

P  =  Pn  +  Pi -h  PplV^  (11) 


where 

Pn  and  Pi  are  the  partial  pressures  of  neutrals  and  ions,  re¬ 
spectively,  in  the  discharge  column; 

Ppi  is  the  mass  plasma  density; 

V  is  the  plasma  velocity. 

In  the  cavity,  the  plasma  velocity  changes  from  zero  near  the 
anode  up  to  a  value  equal  to  the  local  velocity  of  sound  at  the 
cavity  outlet  section.  The  considered  plasma  flow  problem  in¬ 
side  the  cavity  in  principle  is  two-dimensional;  however,  in  the 
main  part  of  the  plasma  bulk  this  axial  velocity  is  small  [26]. 
In  our  pressure  balance,  we  will  use  the  plasma  velocity  as  a 
parameter.  Possible  implications  of  this  assumption  will  be  dis¬ 
cussed  later.  For  known  equilibrium  pressure,  one  can  calculate 
the  ablation  flux  using  the  Langmuir  law  [27] 

where  Tg  is  the  gas  temperature.  Our  estimation  shows  that 
under  considered  conditions  in  the  plasma,  the  bulk  gas  is  in 
thermodynamic  equilibrium  with  the  ions  and  electrons  and  thus 
Tg  =  r,  whereas  near  the  surface  the  gas  has  a  temperature 
equal  to  the  Teflon  surface  temperature.  Thus,  for  calculation 
of  the  ablation  rate,  we  will  use  Tg  —  T^.  Generally,  in  the 
boundary  layer  near  the  Teflon  surface,  the  temperature  sharply 
changes  from  the  plasma  bulk  temperature  to  the  surface  tem¬ 
perature.  The  temperature  profile  in  this  boundary  layer  is  ap¬ 
proximately  linear  [21].  The  characteristic  thickness  of  this  tem¬ 
perature  layer  is  about  an  electron-ion  collision  mean  free  path, 
which  is  about  10“^  m  in  the  considered  range  of  parameters 
and  much  smaller  than  the  cavity  radius. 

The  total  mass  of  the  dielectric  material  ablated  during  the 
entire  discharge  pulse  may  be  calculated  by  integration  of  the 
net  flux  of  ablation  and  deposition  caused  by  back  flux  from  the 


plasma.  In  equilibrium,  the  net  ablated  mass  may  be  determined 
as  follows: 


Ma  =  2TRaL  /  T{T)dt 


I" 


(13) 


where  tp  is  the  pulse  duration. 

For  known  equilibrium  pressure  and  electron  temperature, 
one  can  calculate  the  chemical  plasma  composition  using 
models  developed  previously  [28],  [29].  In  the  considered 
range  of  electron  temperature  (1-2  eV)  and  plasma  density 
(lO^^-lO^"^  ni“^),  we  will  assume  that  polyatomic  molecules 
C2F4  fully  dissociate  [29]  and  we  will  start  our  consideration 
from  the  point  when  we  have  gas  containing  C  and  F.  The 
conservation  of  nuclei  for  the  case  of  Teflon  reads 


2{nc  +  nj)  =  +  (14) 

where  nc  and  nj?  are  the  density  of  neutrals  and  nj  and 
are  the  densities  of  ions,  respectively.  In  LTE,  the  density  of 
electrons  rig,  ions  n"*"  and  neutrals  rin  depend  on  the  plasma 
temperature  according  to  the  Saha  equation  [23] 

=  b(^^  r=*/2  exp(-Jc/r)  =  Ac  (15a) 
\9n/c 

=  T^f^exp{-lF/T)  =  Af  (15b) 

\9n  /  p 

where 

B  is  a  constant; 

Qi  and  ga  are  statistical  weights  of  ion  and  electrons; 

Ic,  If  are  the  ionization  potentials  of  C  and  F,  respec¬ 
tively. 

In  our  case,  the  plasma  consists  of  the  electrons,  neutrals  C  and 
F  and  ions  and  F"*" .  Thus,  we  have  five  unknown  partial  pres¬ 
sures.  Therefore,  an  equilibrium  composition  can  be  calculated 
by  solving  (11),  (14),  and  (15)  supplemented  with  the  quasineu¬ 
trality  condition 


rie 


—  4”  ^P' 


(16) 


From  the  system  of  equations  (13)-(15)  it  is  straightforward  to 
find  an  equation  for  the  electron  density  in  an  explicit  form 

2n^(ne  -f-  Ap)  -h  n^(ne  -j-  Ac)  -  {rih  -  2ne) 

X  {2Ap{nc  +  Ac)  H-  Aci^c  +  ^f)}  =  0  (17) 

where  n/,  is  the  total  heavy  particle  density. 

The  heat  deposited  during  the  discharge  pulse  accumulates 
in  the  skin-layer  6  =  {Xtp/pCp)^-^  and  propagates  within  the 
Teflon  after  the  end  of  the  discharge  pulse  (ip  ^  10  •  ps  in 
PPT-4,  [31]).  During  the  after-pulse  cooling,  the  Teflon  surface 
temperature  can  be  estimated  as  follows  [8],  [30]: 

Z=T,  +  -  To)  (18) 

where  Tgp  is  the  Teflon  surface  temperature  at  the  end  of  the 
discharge  pulse. 


m.  Results 

In  this  section,  we  will  present  results  of  calculation  of  the 
plasma  parameters  during  the  discharge  pulse  in  the  PPT-4 
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Fig.  2.  Plasma  temperature  time  distribution  with  plasma  velocity  as  a  parameter 


developed  recently  at  the  University  of  Illinois  Urbana-Cham- 
paign  [9],  [31].  The  set  of  the  equation  consists  of  nonlinear 
algebraic  equations  and  a  differential  equation  describing  the 
energy  balance.  Equations  (1)-(17)  form  a  set  of  eight  equations 
with  the  following  unknowns:  T,  T, ,  Ua ,  ,  nc ,  tif  j  • 

The  differential  equation  was  integrated  numerically  using 
the  Runge-Kutta  algorithm.  Linearization  of  the  experimental 
current  waveform  causes  some  kinks  in  the  graphics  presented 
below.  In  order  to  simulate  the  energy  input  during  the  pulse, 
the  experimental  current  waveform  has  been  used  [9],  [31].  The 
results  presented  in  this  section  correspond  to  the  following 
PPT-4  design  and  discharge  parameters:  Ra  =  3.15  mm, 
L  =  8.3  mm,  current  peak  of  about  8  kA,  pulse  duration 
is  about  10  //s.  The  plasma  velocity  v  is  used  as  a  param¬ 
eter  of  the  problem  and  is  normalized  by  the  sound  speed 
Cs  =  (2T/mj)°  The  initial  plasma  density  is  used  as  a 
parameter  before  Teflon  heating.  This  density  lies  in  the  range 
of  10^^ -10^^  m“^  and  weakly  affects  the  results  as  shown 
below. 

The  time  evolution  of  the  plasma  temperature  is  shown  in 
Fig.  2.  Initially,  the  plasma  temperature  sharply  increases  and 
peaks  at  about  2-2 .3  e  V,  dependent  on  the  velocity  v .  Toward  the 
pulse  end,  the  plasma  temperature  decreases  to  1-1.5  eV.  In  the 
experiment  [31]  it  was  measured  that  the  peak  electron  temper¬ 
ature  lies  in  the  range  of  2-2.5  eV,  decreases  to  about  1  eV,  and 
varies  slightly  with  axial  distance  from  the  thruster  exit  plane 
(not  shown).  The  model  predictions  and  experimental  data  have 
similar  trend,  which  indicates  that  the  model  prediction  reason¬ 
ably  agrees  with  experimental  data.  It  should  be  noted,  how¬ 
ever,  that  the  data  were  taken  in  the  plasma  plume  at  different 
axial  distances  from  the  thruster  exit  plane  (up  to  0.2  cm).  Thus, 
the  experimental  and  calculation  time  scales  are  different  with 
a  shift  approximately  equal  to  the  plasma  flight  time  from  the 
cavity  to  the  plane  where  data  were  taken. 

The  electron  density  distribution  during  the  pulse  is  plotted  in 
Fig.  3.  It  can  be  seen  that  the  plasma  density  has  a  maximum  of 
about  10^^  m“^  and  then  decreases  toward  the  discharge  end. 


The  plasma  density  peaks  at  about  3  /iS,  which  corresponds  to 
the  peak  in  the  current  pulse.  One  can  also  see  that  the  high 
plasma  density  corresponds  to  the  smaller  velocity,  which  is  a 
result  of  mass  conservation. 

The  temperature  of  the  Teflon  surface  is  shown  in  Fig.  4.  The 
temperature  sharply  increases  during  the  first  2  //s  of  the  dis¬ 
charge  pulse  and  peaks  at  about  650  K,  which  is  above  the  tem¬ 
perature  of  Teflon  decomposition  (^600  K,  [8]).  It  should  be 
noted  that  the  Teflon  surface  temperature  only  slightly  depends 
on  plasma  velocity  v. 

The  anode  sheath  potential  drop  is  displayed  in  Fig.  5  with 
initial  potential  drop  as  a  parameter.  The  initial  value  of 
Uao  is  a  result  of  introducing  the  initial  plasma  density  as  a 
parameter.  It  can  be  seen  that  Ua  is  negative  and  has  a  maximum 
value  of  about  -6T,  regardless  of  the  initial  value  of  Ua .  Toward 
the  pulse  end,  the  absolute  value  of  the  anode  sheath  potential 
drop  decreases  significantly. 

The  chemical  composition  of  the  plasma  is  shown  in  Fig.  6(a) 
for  the  case  v  =  0.3  and  Fig.  6(b)  for  the  case  v  =  0.5. 
One  can  see  that  all  densities  peak  at  about  3  //s.  The  fluo¬ 
rine  (F)  neutral  density  is  larger  than  the  carbon  (C)  neutral 
density,  because  originally  Teflon  has  composition  C2F4.  How¬ 
ever,  the  difference  between  fluorine  ion  density  and  carbon  ion 
density  is  much  smaller  because  carbon  has  a  smaller  ioniza¬ 
tion  energy  [ionization  energy  of  carbon  is  11.3  eV  (at  the  first 
level),  whereas  fluorine  has  a  corresponding  ionization  energy 
of  17.4  eV].  It  can  be  seen  that  the  ionization  degree  is  higher 
in  the  case  of  larger  velocity. 

The  distribution  of  the  input  energy  density  (by  the  Joule  heat) 
and  output  (radiation,  electron  flux  to  the  anode,  and  electron 
and  ion  convection  to  the  Teflon)  are  shown  in  Fig.  7.  It  can  be 
seen  that  the  Joule  heat  density  peaks  at  about  10^^  W/m^  and 
is  mainly  carried  off  by  particle  convection  and  radiation.  The 
energy  input  is  larger  than  output  initially,  for  the  time  <3  /iS 
(plasma  heating),  whereas  for  the  time  >3  ps,  the  energy  input 
is  slightly  smaller  than  the  output  (plasma  cooling).  This  corre¬ 
sponds  to  the  plasma  temperature  behavior,  as  shown  in  Fig.  2. 
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Fig.  3.  Electron  density  distribution  during  the  dischai^e  pulse  with  plasma  velocity  as  a  parameter. 


Fig.  4.  Teflon  surface  temperature  with  plasma  velocity  as  a  parameter. 

The  ablated  mass  approximately  increases  linearly  with 
cavity  length,  as  shown  in  Fig.  8,  where  velocity  is  used  as  a 
parameter.  It  can  be  seen  that  in  the  case  of  smaller  velocity, 
ablated  mass  is  higher.  For  comparison,  some  experimental  data 
[31]  are  also  shown.  One  can  see  that  generally  good  agreement 
between  predicted  and  measured  ablated  mass  is  found. 

It  should  be  noted  that  ablation  mass  was  measured  over  1000 
shots  with  a  pulse  rate  of  about  1.1  Hz  [31].  Thus,  it  is  pos¬ 
sible  that  the  Teflon  surface  temperature  during  the  time  be¬ 
tween  pulses  may  be  higher  than  room  temperature,  which  may 
increase  the  total  mass  ablation.  However,  calculation  of  the  late 
time  ablation  shows  that  it  consumes  no  more  that  1  /ig  per 
pulse. 

rv.  Discussion 

In  this  section,  we  will  discuss  the  implications  of  several  of 
the  assumptions  used  in  the  model  described  above. 


The  basic  and  strongest  assumption  is  considering  the  plasma 
to  be  uniform  in  the  cavity,  which  means  that  the  plasma  prop¬ 
erties  are  described  by  average  plasma  parameters.  In  general, 
however,  the  problem  is  2-D  with  plasma  velocity  and  den¬ 
sity  variation  in  the  radial  and  axial  directions.  Plasma  velocity 
varies  from  zero  near  the  anode  up  to  the  sound  speed  at  the 
cavity  exit  plane,  and  thus,  the  plasma  density  should  decrease 
toward  the  cavity  exit  plane.  However,  if  ablation  is  approxi¬ 
mately  uniform  along  the  cavity,  the  plasma  velocity  increases 
substantially  only  near  the  exit  (^0.8  L,  as  shown  in  [13],  [26]). 
In  the  plasma  bulk,  the  velocity  is  about  (0.3-0.5)C5,  as  was 
used  in  the  model.  It  was  also  found  that  the  plasma  velocity 
does  not  significantly  affect  the  plasma  density  and  tempera¬ 
ture  distributions  during  the  discharge  pulse.  In  spite  of  its  sim¬ 
plicity,  our  approach  is  able  to  predict  reasonable  behavior  of 
the  plasma  temperature  and  ablated  mass,  which  are  in  reason¬ 
able  agreement  with  an  experiment. 


V  —  0.3  and  (b)  v  =  0.5. 
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Fig.  7.  Energy  input  and  output  time  distribution. 


Fig.  8.  Ablated  mass  as  a  function  of  the  cavity  length  with  velocity  as  a  parameter.  Experimental  data  were  taken  from  [31]. 


The  present  approach  used  an  ablation  model  based  on  the 
equilibrium  vapor  pressure  depending  on  the  Teflon  surface 
temperature.  This  approach  is  similar  to  that  used  in  [25]. 
Because  equilibrium  vapor  pressure  is  sensitive  to  the  exact 
value  of  the  surface  temperature,  the  heat  transfer  becomes 
an  important  issue.  In  general,  the  heat  transfer  to  an  ablated 
surface  is  a  complicated  nonlinear  problem  requiring  numerical 
solution  [18],  [25].  In  our  model,  we  have  used  a  temperature 
profile  [17],  which  was  found  to  be  a  fair  approximation  for 
a  reacting  surface.  This  approximation  for  the  temperature 
profile  is  used  for  estimation  of  the  temperature  gradient  at 
the  Teflon  surface.  However,  it  was  found  that  the  temperature 
gradient  on  the  Teflon  surface  has  only  a  small  effect  on  the 
surface  temperature. 


In  order  to  predict  the  chemical  composition  of  the  plasma, 
the  present  model  employed  the  Saha  equation  assuming  local 
thermodynamic  equilibrium  in  the  plasma  column.  However, 
LTE  may  be  established  only  in  relatively  dense  plasmas 
>10^^  m"^  over  a  time  period  larger  than  the  characteristic 
relaxation  time  for  ionization  and  recombination  in  the  plasma 
(^10”^  s).  Calculation  shows  that  during  the  main  part  of  the 
pulse  the  above  requirements  are  fulfilled.  Toward  the  pulse 
end,  however,  the  plasma  density  significantly  decreases  and 
the  relaxation  time  becomes  comparable  to  the  pulse  duration. 
In  this  case,  the  LTE  approach  predicts  a  plasma  ionization 
degree  that  is  much  higher  than  it  should  be  in  the  real  situa¬ 
tion.  Thus,  the  electron  density  calculated  after  about  8  //s  is 
overestimated  and  should  be  considered  only  as  an  upper  limit. 
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Let  US  estimate  the  effect  of  possible  plasma  and  current  com¬ 
pression  by  the  self-magnetic  field  during  the  discharge  pulse. 
At  the  current  peak  (^3  /is),  the  Hall  parameter  (ratio  of  electron 
cyclotron  frequency  to  electron-ion  collision  frequency)  is  less 
than  unity.  Thus,  the  conductivity  is  determined  by  Coulomb 
collisions  and  the  effect  of  the  self-magnetic  field  on  the  cur¬ 
rent  flow  is  small.  An  estimation  shows  that  for  PPT-4  condi¬ 
tions  the  ratio  of  the  gas  kinetic  pressure  to  the  magnetic  pres¬ 
sure  >  1  during  the  discharge  pulse,  except  at 

the  current  peak  point,  where  the  gas  kinetic  pressure  is  slightly 
larger  than  is  the  magnetic  pressure.  However,  the  self-  magnetic 
field  may  have  substantial  effect  on  the  discharge  processes  in 
the  cavity  with  smaller  radius.  Therefore,  further  development 
of  the  pulse  discharge  model  may  be  associated  with  taking  into 
account  the  self-magnetic  field. 


V.  Conclusion 

An  analysis  of  the  physical  processes  in  the  Teflon  cavity  of 
a  PPT  shows  that  the  energy  of  the  plasma  column  is  carried  off 
by  particle  convection  to  the  dielectric  and  by  radiation.  During 
the  pulse,  the  electron  density  peaks  at  about  10^"^  m“^  and  de¬ 
creases  to  10^^  m“^  toward  the  pulse  end,  whereas  the  elec¬ 
tron  temperature  peaks  at  about  2.2  eV  and  decays  to  1.5  eV. 
Teflon  surface  temperature  peaks  at  about  650  K,  which  is  above 
the  temperature  of  decomposition  (^600  K).  Calculated  elec¬ 
tron  temperature  ablated  mass  were  found  to  be  in  good  agree¬ 
ment  with  experimental  measurements.  The  present  model  al¬ 
lows  physical  time-dependent  boundary  conditions  to  be  formu¬ 
lated  for  the  plasma  jet  expansion  problem. 
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Models  are  presented  for  a  Teflon®-fed,  pulsed  plasma  thruster  from  plasma  generation  to  plume  far  field.  A  one¬ 
dimensional  model  using  local  thermodynamic  equilibrium  is  developed  to  describe  the  plasma  generation,  Teflon 
ablation,  and  nozzle  acceleration  processes.  A  computer  code  that  uses  two  different  particle  methods  is  employed 
to  simulate  the  unsteady  plasma  plume  expansion  processes,  including  charge  exchange  collisions.  Results  are 
generated  for  a  pulsed  plasma  thruster  under  development  at  the  University  of  Illinois.  General  features  of  the 
plume  reveal  substantial  differences  in  the  expansion  dynamics  of  the  charged  and  neutral  species.  In  addition, 
there  is  a  noticeable  difference  between  the  behavior  of  the  carbon  and  the  fluorine  ions.  Direct  comparisons 
are  made  between  simulation  results  and  experimentally  measured  data  for  electron  number  density  and  plasma 
potential.  These  comparisons  indicate  strengths  and  weaknesses  of  the  models. 


Introduction 

ULSED  plasma  thrusters  are  receiving  renewed  attention  be¬ 
cause  of  their  ability  to  produce  small  impulse  bits  with 
high  reliability.'  An  accurate  simulation  of  plumes  from  all  elec¬ 
tric  propulsion  devices  is  required  for  the  assessment  of  space¬ 
craft  integration  effects.  For  a  pulsed  plasma  thruster  (PPT)  using 
solid  Teflon®  as  a  propellant,  the  primary  integration  concern  is 
the  deposition  of  highly  condensible  plume  effluent  on  spacecraft 
surfaces. 

The  numerical  investigation  of  PPT  plumes  presents  many  dif¬ 
ficulties.  The  ablation  process  that  leads  to  the  acceleration  of  a 
high-density  plasma  is  not  well  understood.  The  chemical  compo¬ 
sition  of  the  flow  exhausting  from  the  thruster  is  not  known.  The 
plasma  density  is  significantly  higher  than  that  produced  by  grid- 
ded  ion  thrusters  and  Hall  thrusters.  The  neutral  density  at  the  exit 
places  the  flow  in  the  near-continuum  regime.  Finally,  as  a  result  of 
the  pulsed  nature  of  the  device,  all  flow  processes  are  unsteady. 

In  this  study,  the  focus  is  on  a  pulsed  plasma  thruster  called  PPT-4 
that  is  under  development  at  the  University  of  Illinois.^  The  PPT-4 
is  unusual  in  that  the  thrust  force  is  generated  primarily  by  an  elec¬ 
trothermal  mechanism.  In  terms  of  modeling,  this  is  a  convenient 
place  to  start  because  difficulties  with  the  modeling  of  the  electro¬ 
magnetic  effects  can  be  omitted. 

Details  of  the  PPT-4  thruster  are  first  described.  This  paper  reports 
on  two  modeling  efforts  for  this  device.  The  first  is  concerned  with 
modeling  the  plasma  generation,  Teflon  decomposition,  and  nozzle 
acceleration  processes.  A  one-dimensional  model  is  described.  The 
main  goal  of  this  study  is  to  provide  time-varying  boundary  condi¬ 
tions  for  the  subsequent  computation  of  the  PPT  plume.  The  plume 
modeling  based  on  particle  methods  is  then  described.  Results  con¬ 
cerning  the  general  flow  features  of  the  plume  are  presented.  In 
addition,  direct  comparisons  are  made  between  simulation  and  ex- 
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periment  for  electron  number  density  and  plasma  potential.  These 
comparisons  serve  to  partially  validate  the  modeling  and  also  indi¬ 
cate  areas  where  further  development  is  needed. 

University  of  Illinois  PPT-4 

The  thruster  under  consideration  is  being  developed  at  the  Uni¬ 
versity  of  Illinois  and  is  called  the  PPT-4.  A  schematic  diagram  is 
shown  in  Fig.  1 .  The  PPT-4  has  a  coaxial  configuration  with  a  cen¬ 
tral  anode  that  is  5  mm  in  diameter  and  an  annular  cathode  that  is 
43  mm  in  diameter.  The  two  electrodes  are  connected  with  a  30-deg 
half- angle  nozzle.  The  Teflon  is  ablated  from  a  cylindrical  cavity 
with  a  diameter  of  5  mm  and  a  length  of  3.2  mm  that  sits  just  in 
front  of  the  central  electrode.  The  measured  mass  ablation  is  12  ^g 
per  pulse  at  9  J  of  input  energy.  The  pulse  length  is  approximately 
10  //s.  The  overall  specific  impulse  has  been  measured  as  1000  s. 
More  details  can  be  found  in  Ref.  2. 

Modeling  of  Plasma  Generation 

Pulsed  Discharge  Model 

In  the  PPT-4,  plasma  is  generated  in  the  Teflon  cavity  as  a  result  of 
ionization  of  the  Teflon  ablation  products.  The  main  features  of  the 
discharge  in  the  Teflon  cavity  are  analogous  to  ablation  controlled 
arcs.^“^  During  the  pulse  impulse,  the  plasma  is  heated  by  the  joule 
heat  and  is  cooled  by  radiation,  energy  losses  caused  by  particle  con¬ 
vection  to  the  anode  and  dielectric,  and  ionization.  This  interaction 
of  the  plasma  column  with  the  walls  of  the  cavity  heats  the  walls 
rapidly.  The  temperature  of  the  surface  layer  increases  and  reaches 
the  critical  temperature  of  decomposition.  The  principal  mecha¬ 
nism  of  thermal  degradation  of  the  dielectric  is  the  breaking  of  the 
bonds  in  the  chain,  which  consumes  the  main  portion  of  the  energy 
transferred.  As  a  result  of  heating,  dissociation,  and  ionization  of 
the  products  of  decomposition,  the  electron  density  in  the  cavity  in¬ 
creases.  The  pressure  and  temperature  of  the  plasma  in  the  cavity  de¬ 
termine  the  plasma  composition.  In  this  section  we  briefly  describe 
the  model  of  the  pulsed  electrical  discharge  in  the  Teflon  cavity.  A 
more  detailed  description  of  the  model  can  be  found  in  Ref.  6. 

A  simplified  model  of  the  pulsed  discharge  in  the  dielectric  cavity 
is  based  on  the  following  assumptions:  1)  all  plasma  parameter  are 
uniform  within  the  plasma  column;  that  is,  a  diffuse  type  of  dis¬ 
charge  is  assumed;  2)  the  plasma  is  quasineutral;  and  3)  the  plasma 
column  is  in  local  thermod)mamic  equilibrium  (LTE),  and  the  elec¬ 
trons,  ions,  and  neutrals  establish  a  single  equilibrium  temperature. 
As  a  result  of  the  foregoing  considerations  and  assumptions,  the 
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following  energy  balance  equation  in  the  quasi -neutral  plasma  col¬ 
umn  can  be  written: 

=  Qj  -  Qr  -  Qd  -  Qa  d) 

2  at 

where  is  the  electron  density,  T  is  the  plasma  temperature,  Qjh 
the  joule  heat,  is  the  radiative  emission  energy  loss  that  depends 
on  plasma  density  and  temperature,  Qd  is  the  energy  loss  caused  by 
convection  of  electrons  and  ions  to  the  dielectric  that  depends  on  the 
plasma  temperature  and  the  potential  drop  in  the  electrostatic  sheath 
near  the  dielectric,  and  is  the  energy  loss  caused  by  electron 
convection  to  the  anode  that  depends  on  the  anode  sheath  potential 
drop. 

The  anode  current  is  carried  by  electrons  and  controlled  by  the 
potential  drop  in  the  anode  sheath.  The  potential  drop  in  the  sheath 
depends  on  the  plasma  and  current  density  and  may  have  spatial 
variation.  Assuming  that  the  plasma  column  is  uniform,  generally, 
the  sheath  is  slightly  negative  in  order  to  repel  the  excess  of  the 
thermal  electron  current,  so  that  the  electron  current  to  the  anode 
is  equal  to  the  circuit  current.  Similar  physical  reasons  lead  to  the 
appearance  of  the  electrostatic  sheath  near  the  dielectric  surface  to 
ensure  a  zero  net  current.  The  transient  sheath  formation  near  the 
dielectric  surface  is  determined  by  the  dielectric  permittivity.  In  the 
considered  range  of  electron  density  (10^’ m~^),  the  charac¬ 
teristic  charging  time^  is  less  than  10"^°  s,  which  is  much  smaller 
than  the  characteristic  time  of  the  discharge  parameter  changes. 
Therefore  we  use  a  quasi-steady-state  sheath  model.  Furthermore, 
we  assume  that  the  plasma  density  is  so  high  that  the  sheath  thick¬ 
ness  (of  the  order  of  the  Debye  length)  is  much  smaller  than  the 
electron  Larmor  radius  and  thus  the  effect  of  the  self-magnetic  field 
in  the  sheath  can  be  neglected.  As  in  the  case  of  the  anode  sheath, 
the  sheath  potential  drop  should  be  negative  in  order  to  repel  the 
excess  of  the  thermal  electron  current,  so  that  the  electron  current 
is  equal  to  the  ion  current. 

In  LTE,  the  composition  of  the  plasma  is  found  as  a  function  of 
the  state  of  the  plasma  that  is  determined  by  the  pressure  and  tem¬ 
perature.  The  calculation  of  the  composition  of  the  plasma  consists 
of  solving  a  set  of  Saha  equations  supplemented  by  a  pressure  bal¬ 
ance,  conservation  of  nuclei,  and  quasineutrality.  It  is  assumed  that 
the  plasma  pressure  is  equal  to  the  equilibrium  vapor  pressure  based 
on  the  temperature  of  the  Teflon  surface.®  The  Teflon  surface  tem¬ 
perature  is  calculated  by  the  solution  of  the  heat  transfer  equation. 
This  equation  is  subject  to  the  boundary  conditions  at  the  Teflon 
surface  accounting  for  heat  flux  from  the  plasma  by  radiation  and 
particle  convection  and  heat  losses  caused  by  evaporation  of  the 
surface.  In  the  cavity,  the  plasma  velocity  changes  from  zero  near 
the  anode  up  to  a  value  equal  to  the  local  velocity  of  sound  at  the 
cavity  outlet.  The  considered  plasma  flow  problem  inside  the  cavity 
is,  in  principle,  two  dimensional.  However,  in  the  main  part  of  the 
plasma  bulk,  the  two-dimensional  effect  is  small.^  In  the  cavity,  the 


plasma  velocity  varies  from  zero  near  the  anode  up  to  the  sound 
speed  at  the  cavity  exit  plane,  and  thus  the  plasma  density  should 
decrease  toward  the  cavity  exit  plane.  However,  if  ablation  is  ap¬ 
proximately  uniform  along  the  cavity,  the  plasma  velocity  increases 
substantially  only  near  the  exit  (at  approximately  80%  of  the  cavity 
length,  as  shown  in  Ref.  9).  A  quasi-steady  analysis^  shows  that,  in 
the  plasma  bulk,  the  velocity  is  approximately  (0.3-0.5)Cj,  where 
Cs  is  the  sound  speed.  Therefore  the  plasma  velocity  was  used  in 
the  model  as  a  free  parameter.  It  was  found  that  the  plasma  veloc¬ 
ity  does  not  significantly  affect  the  plasma  density  and  temperature 
distributions  during  the  discharge  pulse. 

Plasma  Acceleration  in  the  Conical  Nozzle 

The  electromagnetic  thrust  component  of  the  PPT-4  is  created 
mainly  near  the  thruster  exit  plane,  where  the  azimuthal  self- 
magnetic  field  and  the  radial  component  of  the  discharge  current 
are  maximal.  However,  the  current  pulse  is  essentially  gone  before 
the  main  plasma  cloud  exits  from  the  cavity.  Thus,  the  main  part  of 
the  plasma  generated  in  the  cavity  accelerates  in  the  conical  nozzle 
by  the  gasdynamic  mechanism.  In  this  section  we  describe  a  quasi- 
one-dimensional  model  of  a  continuum,  single-fluid  plasma  flow  in 
the  conical  nozzle  neglecting  magnetic  field  effects. 

We  consider  a  quasi-neutral  plasma  in  which  the  electrons  and 
ions  behave  as  ideal  gases.  The  plasma  temperature  is  assumed  to  be 
constant.  This  is  based  on  the  idea  that  the  effects  on  the  temperature 
of  gasdynamic  expansion  and  plasma  ohmic  heating  will  almost 
cancel.  The  temperatures  obtained  by  using  our  approach  (between 
1  and  3  eV)  are  in  good  general  agreement  with  values  measured 
experimentally  in  the  plume  The  plasma  flow  is  considered  to  be 
sourceless,  and  plasma  losses  to  the  nozzle  wall  and  wall  evaporation 
are  neglected.  Furthermore,  we  consider  the  situation  in  which  all 
properties  vary  only  along  the  axial  direction.  The  plasma  jet  cross 
section  A(z)  is  determined  by  the  nozzle  geometry  as  follows: 

A(z)  =  Ao[l  +  (z//?o)tan9p  (2) 

where  Aq  is  the  initial  cross  section,  Rq  is  initial  jet  radius  equal  to 
the  cavity  radius,  and  9  is  the  cone  half-angle.  Taking  into  account 
these  considerations,  we  find  the  equation  for  the  plasma  velocity 
becomes 

-  Cj  ^  ^  2  2tan9  1 

V  dz  "  Rq  l  +  {z/Ro)tane  ^ 

where  V  is  the  plasma  jet  axial  velocity. 

Time-Dependent  Nozzle  Exit  Conditions 

Here  we  present  results  of  the  calculation  of  the  plasma  parame¬ 
ters  at  the  exit  plane  of  the  PPT-4.  To  simulate  the  energy  input  dur¬ 
ing  the  pulse,  the  experimental  current  waveform  has  been  used.^ 
Linearization  of  the  experimental  current  waveform  causes  some 
perturbations  in  the  results  presented  later.  The  calculations  pre¬ 
sented  in  this  section  correspond  to  the  following  PPT-4  design  and 
discharge  parameters:  cavity  size  Ra  =  3.15  mm,  L  =  8.3  mm,  noz¬ 
zle  half-angle  9  =  30  deg,  current  peak  '^8  kA,  and  pulse  duration 
'-'10/is. 

During  the  discharge  pulse,  plasma  is  generated  in  the  cavity  as 
a  result  of  Teflon  decomposition  and  then  accelerates  in  the  coni¬ 
cal  nozzle.  The  solution  of  the  energy  and  mass  balance  equations 
provides  time  dependencies  of  the  plasma  temperature  and  compo¬ 
sition  in  the  cavity.  During  the  plasma  flow  in  the  nozzle,  plasma 
composition  and  temperature  remain  constant  while  plasma  veloc¬ 
ity  increases  starting  from  the  sound  speed  at  the  cavity  exit  plane. 
In  the  conical  nozzle,  plasma  acceleration  up  to  a  Mach  number  of 
4  is  calculated. 

The  time  evolution  of  the  plasma  component  densities,  tempera¬ 
ture,  and  velocity  ate  shown  in  Figs.  2a-- 2c.  The  chemical  composi¬ 
tion  of  the  plasma  is  shown  in  Fig.  2a.  One  can  see  that  all  densities 
peak  at  ~6  ^s,  which  corresponds  to  the  peak  of  the  discharge  cur¬ 
rent.  The  fluorine  (F)  neutral  density  is  larger  than  the  carbon  (C) 
neutral  density  because  Teflon  has  a  composition  of  C2F4.  However, 
the  difference  between  fluorine  ion  density  and  carbon  ion  density 
is  much  smaller  because  carbon  has  a  smaller  ionization  energy. 
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Fig.  1  Variations  with  time  of  a)  species  number  densities,  b)  electron 
temperature,  and  c)  velocity  at  the  thruster  exit  plane. 

(The  first  ionization  energy  of  carbon  is  11.3  eV,  whereas  fluo¬ 
rine  has  a  corresponding  ionization  energy  of  17.4  eV.)  Initially,  the 
plasma  temperature  sharply  increases  and  peaks  at  '^2.3  eV.  Toward 
the  pulse  end  the  plasma  temperature  decreases  down  to  ~  1.5  eV. 
These  profiles  are  used  as  time-varying  boundary  conditions  for  the 
plasma  plume  simulations. 

Plume  Model 

A  computer  model  for  PPT  plumes  has  been  proposed  by  Yin 
and  Gatsonis.'®  A  hybrid  fluid-particle  approach  was  employed. 


In  the  present  study,  we  employ  a  similar  approach.  Neutrals  (C 
and  F)  and  ions  (C*^  and  F^)  are  modeled  as  particles.  Particle 
collisions  are  computed  by  using  the  direct  simulation  Monte  Carlo 
method  (DSMC).^  *  Both  momentum  exchange  and  charge  exchange 
collisions  are  simulated.  Momentum  exchange  cross  sections  use 
the  model  of  Dalgamo  et  al.,‘^  and  the  collisicm  dynamics  follows 
the  normal  DSMC  procedures  as  described  in  Ref.  1 1 .  The  charge 
exchange  processes  use  the  cross  sections  proposed  by  Sakabe  and 
Izawa.^^  The  computer  code  uses  separate  grids  for  the  collision 
and  plasma  processes.  The  cells  used  for  the  collisions  are  sized 
according  to  the  local  mean  free  path,  as  is  the  usual  case  in  DSMC 
computations.  The  experimental  facility  backpressure  corresponds 
to  a  density  of  ^  10^®  m~^.  This  is  applied  as  an  auxiliary  condition 
with  which  particles  from  the  thruster  can  collide. 

Acceleration  of  the  charged  particles  in  self-consistertf  electric 
fields  is  simulated  by  using  the  particle-in-cell  method  (PIC).^'*  The 
plasma  potential  0  is  obtained  by  assuming  charge  neutrality  to 
determine  the  electron  number  density  from  the  total  ion  density.  The 
electron  number  density  n,  is  then  used  in  the  Boltzmann  relation 
to  obtain  the  plasma  potential: 

<f>  =  TMnefnrAt)]  (4) 

where  the  electron  temperature  Te  is  in  electron  volts  and  ritef  is  the 
electron  number  density  at  a  reference  point.  This  approach  was 
used  in  our  previous  work  on  Hall  thruster  plumes.*^  In  the  case 
of  the  PPT,  the  reference  point  for  the  Boltzmann  relation  is  taken 
as  the  thruster  exit.  It  is  assumed  that  the  potential  here  is  constant. 
The  variation  of  electron  number  density  at  the  thruster  exit  obtained 
in  the  plasma  generation  modeling  is  used  to  change  the  reference 
density  as  a  function  of  time.  A  constant  electrcm  temperature  of 
2  eV  is  used  in  the  Boltzmann  relation  throughout  the  computation. 
Because  charge  neutrality  is  assumed,  the  PIC  cells  do  not  have  to 
be  of  the  order  of  the  Debye  length.  Instead  they  are  chosen  to  be 
small  enough  to  resolve  in  a  reasonable  way  the  gradients  in  the 
potential.  The  grids  used  in  the  computation  are  shown  in  Fig.  3. 
A  single  time  step  given  by  the  reciprocal  of  the  maximum  plasma 
frequency  is  used  throughout.  All  results  are  time  dependent  and 
are  integrated  over  small  intervals  of  time,  which  are  of  the  order  of 
2  X  10-“^  s. 

The  main  difficulties  in  performing  the  computations  of  the  PPT 
plume  arise  from  the  transient  nature  of  the  expansion.  Because  the 
flow  conditions  change  as  a  function  of  time,  the  characteristic  time 
and  length  scales  of  both  collision  and  plasma  phenomena  continu¬ 
ally  change  throughout  the  computational  domain.  In  principle,  the 
simulations  should  use  grids  that  change  with  the  local,  transient 
flow  conditions.  In  this  first  study,  we  have  adopted  the  simpler  ap¬ 
proach  of  employing  grids  that  should  capture  most  of  the  physics 
most  of  the  time.  Another  numerical  problem  related  to  the  transient 
flow  conditions  is  the  fact  that  there  are  often  regions  of  the  flow 
where  the  number  of  simulated  particles  is  very  small.  This  occurs 
at  the  leading  edge  of  the  forward  expansion  at  early  times  and  at 
the  trailing  edge  of  the  expansion  at  late  times.  The  low  number  of 


Z(m) 

Fig.  3  Computer  grids  employed  in  the  plume  particle  simulation. 
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particles  in  these  regions  leads  to  significant  scatter  in  the  simula¬ 
tion  results  in  these  regions.  In  the  early  stages  of  this  investigation 
we  experimented  with  the  use  of  very  large  computations  involving 
many  millions  of  particles  perform^  on  parallel  computers.  This 
approach  only  moved  the  statistical  scatter  to  lower  density  portions 
of  the  flowfield.  As  shown  in  the  Results  section,  the  density  varies 
by  several  orders  of  magnitude  across  the  jet  expansions  and  so 
the  statistical  problem  is  unavoidable.  It  should  be  noted,  however, 
that  most  of  the  results  presented  are  unaffected  by  the  statistical 
fluctuations. 

Results 

Particles  are  introduced  into  the  flowfield  through  the  nozzle  exit 
plane  by  using  the  output  from  the  plasma  generation  model  (Fig.  2). 
A  uniform  radial  profile  of  properties  is  assumed  across  the  exit 
plane  except  that  a  conical  distribution  of  flow  angle  is  used.  Ions 
and  neutrals  are  assumed  to  have  the  same  velocity  and  the  same 
temperature.  This  assumption  is  based  on  the  fact  that  the  number 
densities  are  relatively  high  so  that  the  frequency  of  charge  exchange 
collisions  inside  the  nozzle  is  high.  Simulations  are  performed  with 
and  without  the  facility  backpressure.  It  is  found  that  the  results  are 
almost  identical,  indicating  that,  unlike  some  other  electric  propul¬ 
sion  devices  such  as  Hall  thrusters  and  ion  thrusters,  the  PPT  plume 
is  not  affected  by  the  chamber  background  gas. 

To  illustrate  some  of  the  basic  dynamics  of  the  time  evolution  of 
the  PPT  plume  expansion  processes.  Figs.  4a-4d  show  contours  of 
plasma  potential  at  four  different  times.  In  these  results  it  is  clear 
that  the  plasma  potential  expands  very  rapidly  about  a  source  that 
moves  with  the  plume.  The  gradients  in  potential  are  such  as  to  create 


electric  fields  that  will  accelerate  ions:  1)  in  the  forward  direction  at 
the  front  of  the  expansion;  2)  in  the  backward  direction  at  the  rear 
end  of  the  expansion,  which  results  in  backflow  of  ions  behind  the 
thruster,  and  3)  in  the  radial  direction  away  from  the  axis.  Of  course, 
these  forces  act  only  on  the  charged  particles;  thus  in  the  results  that 
follow,  it  is  expected  that  the  ions  will  behave  quite  differently  from 
the  neutral  atoms. 

To  demonstrate  the  different  dynamics  of  ions  and  neutrals. 
Figs.  5a-5d  show  the  variation  of  all  species  densities  along  the 
axis  at  four  dilferent  times.  As  shown  in  Fig.  5a,  at  10  ^s  after 
ignition,  the  species  densities  follow  the  profiles  shown  in  Fig.  2a. 
The  ions  are  in  a  larger  concentration  than  the  neutral  atoms,  and 
fluorine  is  more  abundant  than  carbon.  The  relative  concentrations 
of  the  species  do  not  change  dramatically  at  later  times  as  shown 
in  Figs.  5b-5d.  However,  the  shape  of  the  profiles  reveals  that  the 
charged  species  are  spread  over  a  larger  axial  region  than  the  neu¬ 
trals.  The  neutral  atoms  expand  as  a  relatively  compact  puff  of  gas, 
whereas  the  ions  are  spread  out  both  ahead  of  and  behind  the  center 
of  expansion  as  a  result  of  electric  field  effects. 

Comparisons  between  the  computation  and  experimental 
measurement^  of  the  electron  number  density  along  the  axis  are 
shown  at  three  different  locations  in  Figs.  6a-6c.  Because  charge 
neutrality  is  assumed,  the  computational  results  are  obtained  by 
summing  the  carbon  and  fluorine  ion  densities.  When  a  compari¬ 
son  is  made  with  the  experimental  data,  it  is  expected  that  the  data 
collected  over  the  first  10  //s  after  ignition  cannot  be  used.  This  is 
partly  because  of  the  noise  introduced  into  the  measurements  by  the 
operation  of  the  spark  plug  igniter.  In  addition,  there  is  no  attempt  in 
our  modeling  to  simulate  the  behavior  of  the  ignition  plasma.  With 


a)  10,  b)  15,  c)  20,  and  d)  25  0  s  after  ignition. 
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Fig.  5  Profiles  of  species  number  density  along  the  axis;  a)  10,  b)  15,  c)  20,  and  d)  25  0  s  after  ignition. 


this  in  mind,  we  can  make  two  clear  observations  about  the  com¬ 
parisons  shown  in  Fig.  6.  The  first  is  that  there  is  a  tendency  in  the 
simulations  to  overpredict  the  density  peak.  However,  it  should  be 
noted  that  the  experimental  uncertainty  is  ±50%  and  the  predicted 
peaks  lie  within  this  range.  There  also  appears  to  be  some  indication 
of  saturation  in  the  experimental  measurements.  It  is  surprising  to 
see  such  a  small  reduction  in  the  peak  density  measured  at  10  and 
14  cm  from  the  thruster.  The  second  clear  conclusion  is  that  the 
simulation  underpredicts  the  densities  at  both  early  and  late  times. 

A  comparison  of  the  simulation  results  with  experimental  data 
shows  that  the  main  plasma  plume  features  are  captured,  although 
the  simulation  underestimates  the  densities  at  long  times.  The  un¬ 
derprediction  of  the  plasma  density  at  long  times  is  a  result  of  the 
basic  assumption  of  the  plasma  generation  model  in  considering 
the  plasma  to  be  uniform  in  the  cavity.  This  approach  means  that 
the  plasma  properties  are  described  by  average  plasma  parameters. 
Thus,  this  model  predicts  that  the  plasma  will  leave  the  cavity  after 
the  pulse  end  in  a  time  that  is  short  compared  with  the  pulse  duration. 
In  general,  however,  the  problem  is  two  dimensional,  with  plasma 
velocity  and  density  variation  in  the  radial  and  axial  directions.  The 
plasma  velocity  varies  from  zero  near  the  central  electrode  up  to  the 


sound  speed  at  the  cavity  exit  plane.  After  the  pulse  ends,  there  is  a 
substantial  reduction  in  the  ablation  rate  and  the  plasma  density  in 
the  cavity  will  decrease.  Plasma  density  temporal  decay  in  this  case 
will  be  much  slower  than  that  in  the  present  model.  In  the  future,  this 
effect  will  be  included  in  our  model  by  considering  spatial  variation 
in  the  cavity  properties. 

To  predict  the  chemical  composition  of  the  plasma,  the  present 
model  uses  the  Saha  equation  assuming  local  thermodynamic  equi¬ 
librium  in  the  plasma  column.  However,  LTE  may  be  established 
only  in  relatively  dense  plasmas  with  an  electron  number  density 
greater  than  10^“  m”^  for  a  time  period  larger  than  the  characteris¬ 
tic  relaxation  time  for  ionization  and  recombination  in  the  plasma 
(of  the  order  of  lO"”^  s;  Ref.  16).  Calculations  show  that  during  the 
main  part  of  the  pulse,  these  requirements  are  fulfilled.  Toward  the 
pulse  end,  however,  the  plasma  density  significantly  decreases  and 
the  relaxation  time  becomes  comparable  with  the  pulse  duration. 
Thus,  toward  the  pulse  end,  the  LTE  approach  predicts  a  plasma 
ionization  degree  that  is  much  higher  than  it  should  be  in  the  real 
situation. 

In  Figs.  7a-7c  a  comparison  is  made  between  measurements  and 
computations  for  the  plasma  potential  on  the  axis  at  the  same  three 
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locations  as  considered  in  Fig.  6.  If  we  again  ignore  the  first  10 
we  can  see  that  the  initial  rise  and  the  actual  peak  of  the  poten¬ 
tial  data  at  all  three  locations  is  well  predicted  by  the  simulation. 
Also,  consistent  with  the  electron  number  density  comparisons,  the 
simulations  show  a  more  rapid  decay  in  potential  at  long  times  after 
ignition.  Overall,  the  comparisons  with  the  experimental  data  shown 
in  Figs.  6  and  7  indicate  that  the  modeling  has  achieved  a  reasonable 
level  of  success. 

As  mentioned  earlier,  the  primary  reason  for  performing  plume 
computations  of  the  PPT  is  to  assess  possible  spacecraft  interaction 
effects.  With  this  in  mind,  we  present  in  Figs.  8  and  9  mass  flux  re¬ 
sults  as  a  function  of  time  from  the  simulations  in  the  radial  planes 
at  50  cm  forward  of  the  thruster  exit  and  above  the  thruster  exit, 
respectively.  The  latter  plane  is  chosen  to  assess  the  potential  back- 
flow  problem.  The  forward  contamination  problem  could  occur  if 
PPTs  are  employed  on  closely  spaced  spacecraft  flying  in  formation. 
Figure  8  indicates  that  the  mass  flux  at  any  plane  changes  dramat¬ 
ically  as  a  function  of  time.  As  shown  in  Fig.  8a,  the  earliest  flux 
consists  of  the  fastest  ions  that  are  accelerated  electrostatically  ahead 
of  the  neutrals.  Because  the  charge-to-mass  ratio  of  carbon  is  about 
50%  larger  than  that  for  fluorine,  it  is  the  carbon  ions  that  are  the 


first  species  to  impact  a  surface  forward  of  the  thruster.  Figure  8b 
indicates  that  any  significant  flux  of  neutral  atoms  is  delayed  by 
~5  ^s  in  reaching  the  surface  and  even  then  the  flux  magnitude  is 
~  2  orders  of  magnitude  below  that  of  the  ions.  At  later  times,  the 
radial  acceleration  of  ions  leads  to  their  relative  magnitudes  of  mass 
flux  becoming  comparable  with  that  of  fluorine  atoms.  The  mass 
flux  of  carbon  atoms  remains  at  a  lower  level  up  until  35  fis. 

The  backflow  contamination  issue  is  demonstrated  in  Fig.  9, 
and  this  phenomenon  would  be  of  concern  on  any  spacecraft  us¬ 
ing  PPTs.  Relative  to  ion  thrusters  and  Hall  thrusters,  the  poten¬ 
tial  for  backflow  with  PPTs  is  significantly  higher.  This  occurs 
because  the  carbon  and  fluorine  ions  are  significantly  more  mo¬ 
bile  than  xenon  and  because  the  neutral  and  ion  densities  are  orders 
of  magnitude  higher  for  the  PPT,  leading  to  increased  collisional 
scattering.  In  the  simulations,  almost  no  backflow  of  neutral  atoms 
is  predicted.  The  peak  mass  flux  occurs  soon  after  ignition  as  shown 
in  Fig.  9a,  and  the  maximum  level  is  only  1  order  of  magnitude 
below  the  maximum  found  in  Fig.  8  for  the  forward  fluxes.  Again 
because  of  their  increased  mobility,  carbon  ions  dominate  the  back- 
flow  flux.  At  later  times,  it  can  be  seen  that  the  mass  flux  decreases 
dramatically. 
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Fig.  9  Backflow  mass  flux  in  the  plane  above  the  thruster  exit:  a)  10,  b)  15,  c)  20,  and  d)  25  ^  s  after  ignition. 


Conclusions 

A  model  has  been  developed  to  describe  the  plasma  processes 
of  a  Teflon-fed  pulsed  plasma  thruster  from  plasma  generation  to 
plume  far  field.  The  device  considered  was  the  PPT-4,  which  is  being 
developed  at  the  University  of  Illinois.  From  the  modeling  stand¬ 
point,  this  is  an  interesting  thruster,  as  the  acceleration  is  generated 
primarily  by  electrothermal  effects. 

The  plasma  generation,  Teflon  ablation,  and  nozzle  expansion 
were  modeled  by  using  a  one -dimensional  approach.  Local  thermo¬ 
dynamic  equilibrium  was  used  to  compute  the  chemical  composi¬ 
tion.  A  quasi-steacly  assumption  was  made  to  compute  the  nozzle 
flow.  The  time-dependent  properties  computed  at  the  nozzle  exit 
were  used  as  boundary  conditions  for  a  separate  plume  computation. 

The  plume  computation  used  a  combination  of  the  PIC  and  DSMC 
method  to  compute  the  plasma  and  collision  phenomena  in  the  ex¬ 
panding  plume.  The  computations  of  plasma  potential  indicated  that 
the  ions  are  accelerated  ahead  of  the  expanding  plume,  radially  away 
from  the  axis,  and  backward  behind  the  thruster.  It  was  found  that, 
as  a  result  of  electric  field  effects,  the  dynamics  of  the  ions  and  that 
of  the  atoms  in  the  plume  were  substantially  different.  In  addition, 
because  of  a  relatively  large  difference  in  the  charge -to-mass  ra¬ 
tio,  the  dynamics  of  the  carbon  and  that  of  the  fluorine  ions  were 


also  perceptibly  different.  It  was  found  that  the  primary  source  of 
backflow  contamination  from  the  thruster  was  carbon  ions.  The  for¬ 
ward  contamination  was  found  to  be  time  dependent,  with  an  initial 
flux  of  ions  followed  by  a  later  flux  of  neutral  atoms. 

Comparisons  were  made  between  the  computations  and  exper¬ 
imental  measurements  taken  at  the  University  of  Illinois  for  the 
plasma  potential  and  electron  number  density.  In  general,  the  agree¬ 
ment  between  the  data  sets  was  very  good.  It  was  found  that,  at 
late  times,  the  computations  predicted  a  more  rapid  decay  in  both 
potential  and  electron  number  density  in  comparison  with  the  mea¬ 
sured  data.  This  suggested  improvements  in  the  plasma  generation 
model  that  will  be  investigated  further.  The  good  agreement  ob¬ 
tained  between  predictions  and  measurements  provides  confidence 
in  the  overall  modeling  approach  used  in  this  study. 
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Abstract 

Progress  in  combined  device/plume  modeling  is  presented  for  a  Teflon-fed,  pulsed  plasma 
thruster  from  plasma  generation  to  the  plume  far  field.  In  this  work  we  apply  a  one-dimensional 
unsteady  model  for  the  plasma  generation  and  acceleration  process.  A  new  kinetic  ablation 
algorithm  is  employed  to  calculate  the  Teflon  ablation  rate  as  a  function  of  plasma  parameters.  A 
near  cathode  sheath  model  is  included  to  calculate  the  plasma  potential  at  the  thruster  exit  plane. 
Results  are  compared  with  data  for  the  electrothermal  device,  PPT-4.  Performance 
characteristics  of  the  PPT  such  as  mass  ablation  and  thrust  impulse  are  calculated.  Predicted 
plasma  properties,  thruster  performance  and  plasma  parameter  distribution  in  the  plume  are 
found  to  be  In  agreement  with  available  experimental  data. 

Introduction 

Pulsed  plasma  thrusters  (PPT’s)  have  combined  advantages  of  system  simplicity,  high  reliability, 
low  average  electric  power  requirement  and  high  specific  impulse'.  The  PPT  is  considered  as  an 
attractive  propulsion  option  for  orbit  insertion,  drag  makeup  and  attitude  control  of  small  satellites. 
PPT’s,  however,  have  very  poor  performance  characteristics  and  an  overall  efficiency  at  the  level 
of  about  To  improve  the  PPT  performance  several  directions  are  being  considered^ 

Accurate  simulation  of  these  devices  and  plumes  is  required  for  the  design  of  PPT’s  with 
improved  performances  and  for  assessment  of  spacecraft  integration  effects. 
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In  the  present  study  we  concentrate  on  the  pulsed  plasnna  thruster  called  PPT-4  that  was 
developed  recently  at  the  University  of  liiinois^  This  is  an  electrothermal  device  that  derives  most 
of  its  acceleration  from  the  eiectrothermal  or  gasdynamic  mechanism.  This  thruster  is  axially 
symmetric  and  a  discharge  occurs  between  the  annular  cathode  at  the  thruster  exit  plane  and  the 
circular  anode  located  at  the  far  end  of  a  cylindrical  cavity  made  of  Teflon.  The  plasma  generated 
inside  this  cavity  is  acceierated  in  a  diverging  nozzle  that  is  attached  to  the  downstream  end  of 
the  cavity.  The  device  has  a  pulse  length  of  about  10  ps.  and  the  overall  specific  impulse  was 

measured  to  be  850  s. 

in  a  series  of  previous  papers,  we  describe  our  efforts  to  model  various  aspects  of  this 
electrothermal  PPT.  In  Ref.  5,  a  model  of  the  Teflon  ablation  and  plasma  discharge  processes  is 
described.  The  model  was  calibrated  against  mass  ablation  data  from  the  PPT-4.  In  Ref.6,  the 
charging,  heat  and  flow  effects  associated  with  large  Teflon  particulates  in  the  plasma  jet  of  the 
electrothermal  PPT  were  considered.  Here  it  was  predicted  that  the  small  macro-particles  are 
expected  to  decompose  within  the  plasma  jet.  Finally,  in  Ref.  7,  the  results  obtained  in  Ref.  5,  at 
the  thruster  nozzle  exit  were  used  as  boundary  conditions  to  perform  a  particle-based  PIC-DSMC 
computation  of  the  electrothermal  PPT  plume.  A  significant  conclusion  from  Ref.  7  was  that 
almost  all  of  the  back-flow  to  the  spacecraft  from  this  device  arises  from  carbon  ions  due  to  their 

high  mobility. 

The  main  physical  processes  in  this  type  of  PPT  occur  in  the  Teflon  cavity.  Rapid  heating  of  a  thin 
dielectric  surface  layer  leads  to  decomposition  of  the  material  of  the  wall.  As  a  result  of  heating, 
decomposition  and  partial  ionization  of  the  decomposition  products,  the  total  number  of  particles 
increases  in  the  cavity.  The  problem  of  the  ablated  controlled  discharge  has  a  generai  interest 
since  it  can  be  used  for  various  appiications®-®  '°.  In  these  devices,  the  discharge  energy  is 
principally  dissipated  by  ablation  of  wall  material,  which  then  forms  the  main  component  of  the 
discharge  plasma.  The  ablated  vapor  increases  the  pressure  within  the  capiliary  and  the  plasma 
is  expelled  through  the  exit.  Previously,  discharge  evolution  in  the  PPT-4  Tefion  cavity  was 
studied  by  Keidar  et  aP  assuming  uniform  plasma  parameters.  However,  further  understanding 
of  the  physical  processes  involved  requires  more  detailed  analyses  including  spatial  variation  of 
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the  plasma  parameter  distribution  along  the  cavity  and  a  more  sophisticated  ablation  model.  The 
present  model  of  the  capillary  discharge  employs  a  recently  developed  ablation  model'^  These 
changes  allow  us  to  provide  more  accurate  boundary  conditions  for  the  plume  simulation. 

The  capillary  discharge  model 

The  model  presented  here  describes  the  physics  of  the  plasma  generation  and  acceleration  in  a 
Teflon  cavity  for  a  pulsed  electrical  discharge  as  shown  in  Fig.  1.  The  mam  features  of  the  model 
of  the  electrical  discharge  in  the  dielectric  cavity  include  Joule  heating  of  the  plasma,  heat 
transfer  to  the  dielectric.  Teflon  ablation  and  electrothermal  acceleration  of  the  plasma  up  to  the 
sound  speed  at  the  cavity  exit.  Mechanisms  of  energy  transfer  from  the  plasma  column  to  the  wall 
of  the  Teflon  cavity  includes  heat  transfer  by  particle  convection  and  by  radiation.  The  Teflon 
ablation  is  based  on  a  recently  developed  kinetic  ablation  model’^  It  is  assumed  that  all 
parameters  vary  in  the  axial  direction  x  (see  Fig.  1).  Since  the  axial  pressure  and  velocity 
gradients  are  much  greater  than  the  radial  gradients  we  assume  that  radial  variation  of  plasma 
temperature,  pressure  and  velocity  are  neglegible’^’^.  The  axial  component  of  the  mass  and 
momentum  conservation  equations  read: 

ap/at  +  a(pv)  =  . 

av/at +vav/ax  = . . 


where  p  is  the  plasma  density,  P  is  the  pressure,  V  is  the  plasma  velocity  and  r(t,x)  is  the 
ablation  rate.  The  energy  balance  equation  can  be  written  in  the  form: 


|ne  (aT/at  +  VzaT/az)  =  Qj  -  Qr  -  . . . . 

where  Qj  is  the  Joule  heat,  Qr  is  the  radiation  heat  and  Qf  is  the  heat  associated  with  particle 
fluxes.  This  equation  depends  on  the  coordinate  along  the  cavity.  However,  our  estimation  and 


3 


previous  calculations  show’"  that  the  arc  temperature  varies  only  slightly  with  axial  position  and 
therefore  we  further  assume  aT/9z=0.  The  Teflon  surface  temperature  is  calculated  from  the  heat 
transfer  equation  with  boundary  conditions  that  take  into  account  vaporization  heat  and 
conductivity.  The  solution  of  this  equation  is  considered  for  two  limiting  cases  of  substantial  and 
small  ablation  rate  very  similar,  to  that  described  in  Ref.  11.  For  known  pressure  and  electron 
temperature  one  can  calculate  the  chemical  plasma  composition  assuming  LTE”’’®'’®.  The  Saha 
equations  are  supplemented  by  the  conservation  of  nuclei  and  quasi-neutrality. 


Electrostatic  sheaths 

The  electrostatic  sheath  near  the  cathode  provides  the  current  continuity  from  the  cathode  to  the 
plasma  bulk  as  shown  in  Fig.  1.  We  assume  that  the  cathode  emission  plays  a  small  role  in  the 
current  balance.  The  total  current  density  in  the  sheath  consists  of  electron  Je  and  ion  J,  current 

densities; 

. (4) 

J  =  Jj  + Je . 

In  the  case  of  a  planar  sheath  in  front  of  the  cathode  (the  Debye  radius  Is  much  less  than  the 
cathode  length)  Jj  is  determined  by  the  Bohm  relation: 

. (5) 

Jj  =  . . 

where  n  is  the  plasma  density  at  the  plasma-sheath  interface  (see  Fig.  1).  The  electron  current  is 
due  to  high  energy  electrons  that  penetrate  the  electrostatic  barrier: 

/0\ 

jg  =  %  . . 

where  M  is  the  potential  drop  across  the  near-cathode  sheath.  For  given  current  density  one  can 
calculate  the  potential  drop: 

Act)  =  Te  ln((me/mi)°^  -  . . 

One  can  see  that  the  potential  drop  depends  upon  current  density,  plasma  density  and  electron 
temperature. 
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In  the  cavity  near  the  Teflon  surface,  the  electrostatic  sheath  potential  drop  is  negative  in  order  to 
repel  the  excess  thermal  electrons,  so  that  the  electron  current  Je  is  equal  to  the  ion  current  J,  It 
was  concluded  previously  that  during  the  discharge  pulse,  a  quasi-steady  sheath  structure  is 
formed  and  that  under  typical  PPT  conditions  this  sheath  is  unmagnetized  in  the  self-magnetic 
field  generated  during  the  pulse’’.  Under  the  conditions  mentioned  above,  the  potential  drop  in 

the  sheath  is  calculated  as: 


Ud  =  -T  In  (Jelh^Ji) 


(8) 


Where  Je.h  is  the  random  electron  current  density  and  Ji  is  the  ion  current  density  also  determined 


by  the  Bohm  condition. 


Ablation  model 


The  ablation  model  employed  here  is  based  on  a  kinetic  model  of  the  Knudsen  layer  near  the 
ablated  surface,  which  was  analyzed  using  the  distribution  function  moment  method  .  This 
method  employs  an  approximation  of  the  distribution  function  within  the  non-equilibrium  Knudsen 
layer  as  a  sum  of  the  distribution  functions  before  and  after  this  layer  with  a  coordinate  dependent 
coefficient.  In  our  problem  of  evaporation,  it  is  only  important  to  know  the  parameters  on  the 
boundaries  and  not  their  variation  between  the  boundaries.  This  means  that  the  problem  is 
reduced  to  the  integration  of  the  conservation  equations  of  mass,  momentum  and  energy: 


JVxf(V)  dV  =  const 

JVx^f(V)dV  =  const . 

jVxV^  f(V)dV  =  const 

After  integration  of  equation  (9)  we  obtain  a  set  of  equations  in  which  parameters  at  the  external 
boundary  of  the  Knudsen  layer  depend  upon  velocity  at  the  Knudsen  layer  edge.  Applying  mass 
and  momentum  conservation  between  the  edges  of  the  hydrodynamic  layer,  one  can  find  the 
velocity  at  the  outer  boundary  at  the  Knudsen  layer.  Velocity  and  density  at  the  outer  boundary  of 
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the  Knudsen  layer  determine  the  ablation  rate.  The  system  of  equations  Is  closed  If  Ihe 
equilibrium  vapor  pressure  can  be  specified.  In  Ihe  case  of  Teflon,  the  equilibrium  pressure 

formula  is  used’; 

P  =  . . 

where  P  is  the  equilibrium  pressure.  Pc  and  Tc  are  the  characteristic  pressure  and  temperature. 

respectively. 


Nozzle  and  Plume  models 

The  plasma  flow  through  the  conical  nozzle  is  modeled  by  a  quasi-one-dimensional  continuum 
approach  similar  to  that  used  previously  (Ref.  7).  We  have  assumed  that  the  main  part  of  the 
plasma  generated  in  the  cavity  accelerates  in  the  nozzle  by  the  gasdynamic  mechanism.  We 
have  considered  the  quasi-neutral  plasma  where  ions  and  electrons  are  assumed  to  be  ideal 
gases.  This  model  also  relies  on  the  assumption  that  the  flow  is  sourceless  and  plasma  losses  to 
the  wall  and  wall  evaporation  are  small  and  can  be  neglected. 

The  plume  model  is  based  on  a  hybrid  fluid-particle  approach  similar  to  that  used  previously 
(Refs.  7.  20).  In  this  model,  the  neutrals  and  ions  are  modeled  as  particles  while  electrons  are 
treated  as  a  fluid.  Elastic  (momentum  transfer)  and  non-elastic  (charge  exchange)  collisions  are 
included  in  the  model.  The  particle  collisions  are  calculated  using  the  direct  simulation  Monte 
Carlo  (DSMC)  method^’.  Momentum  exchange  cross  sections  use  the  model  of  Dalgarno  et  al.  , 
while  charge  exchange  processes  use  the  cross  sections  proposed  by  Sakabe  and  Izawa"^ 
Acceleration  of  the  charged  particles  in  the  self-consistent  electric  fields  is  computed  using 
Particle-ln-Cell  method  (PIC)^'*.  We  have  assumed  quasi-neutrality  that  allows  determination  of 
the  electron  density.  The  plasma  potential  with  respect  to  the  thruster  exit  plane  is  calculated 
using  the  Boltzmann  relation.  The  plasma  potential  at  the  thruster  exit  plane  with  respect  to  the 
cathode  varies  with  time  and  is  calculated  using  the  near-cathode  sheath  model.  The  grids 
employed  in  this  computation  are  similar  to  those  used  previously  (Ref. 7). 
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Results 


In  this  section  we  present  results  of  the  calculation  of  the  plasma  parameter  temporary  and 
spatial  variation  in  the  Teflon  cavity  of  the  PPT-4  and  also  plume  flowfield  based  on  the  boundary 
conditions  provided  by  the  device  model.  In  addition,  some  thruster  performance  characteristics 
(thrust  impulse  and  ablation  mass)  are  calculated  and  compared  with  available  experimental 
data.  In  this  paper  we  also  show  how  the  improved  thruster  model  and  boundary  conditions  at  the 
thruster  exit  plane  affect  the  plume  simulation. 

Firstly,  we  present  results  for  the  thruster  device  modeling.  These  results  later  will  be  employed 
as  boundary  conditions  for  the  plume  simulation.  All  results  are  presented  for  PPT-4  with  the 
following  geometry  and  discharge  parameters;  anode  radius  Ra=3  mm  and  cavity  length  L=8.3 
mm,  current  peak  of  about  8  kA  and  pulse  duration  of  about  10  ps  (Refs.  4, 25). 

To  illustrate  the  effect  of  the  new  ablation  model  on  the  ablation  rate  calculation,  the  trajectory  of 
the  ablation  rate  in  the  plasma  density  -  Teflon  surface  temperature  plane  during  the  PPT-4  pulse 
is  shown  in  Fig.  2.  One  can  see  that  the  ablation  rate  peaks  at  about  110  kg/m^s  at  3  ps  and  then 
rapidly  decreases.  It  should  be  noted  that  in  the  experiment  the  average  ablation  rate  was 
estimated  to  be  about  30  kg/m^s  (Ref.  4). 

The  time  evolution  of  the  temperature  and  ionization  degree  (at  the  exit  plane,  x=L)  is  shown  in 
Fig.  3.  It  can  be  seen  that  the  electron  temperature  initially  increases  rapidly  and  peaks  at  about  3 
eV  and  then  decreases  to  1  eV  toward  the  pulse  end.  The  model  predicts  that  initially  plasma  in 
the  cavity  is  strongly  ionized  while  after  about  3  ps  the  ionization  degree  decays  substantially. 

To  demonstrate  the  different  ion  and  neutral  species  temporai  and  spatial  variation,  their 
distributions  are  shown  in  Fig.  4.  Both  ion  species  peak  at  early  times  (about  3ps)  while  neutral 
species  peak  later.  The  relative  concentration  of  the  species  changes  along  the  cavity  length  and 
also  during  the  discharge  pulse. 
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The  spatial  and  temporal  distribution  of  the  Teflon  surface  temperature  is  shown  in  Fig.  5.  The 
temperature  sharply  increases  during  the  first  2  ps  of  the  discharge  pulse  and  peaks  at  about  640 
K.  One  can  see  that  the  temperature  varies  only  slightly  along  the  cavity. 

To  illustrate  the  ability  of  the  model  to  predict  some  thruster  performance  characteristics,  we 
calculate  the  thrust  impulse  bit  and  ablation  mass  per  pulse  as  a  function  of  cavity  geometry  and 
compare  with  experimental  data.  The  gasdynamic  thrust  impulse  is  generated  due  to  the 
pressure  force  on  the  anode.  We  integrate  the  pressure  during  the  discharge  to  calculate  the 
thrust  impulse  bit.  In  the  case  of  PPT-4  nozzle  with  the  area  ratio  of  about  100,  the  nozzle  may 
increase  the  thrust  by  a  factor  of  1.5  (Ref.  26).  The  thrust  impulse  bit  dependence  on  the  cavity 
length  is  shown  in  Fig.  6.  One  can  see  that  it  increases  by  a  factor  of  2.5  when  L  increases  from  3 
mm  up  to  25  mm.  A  similar  trend  was  also  found  in  the  experiment^.  The  calculated  Teflon  mass 
ablated  per  pulse  is  shown  in  Fig.  7  as  a  function  of  cavity  length  and  radius.  Generally,  the 
ablated  mass  increases  with  increasing  cavity  length  and  decreasing  cavity  radius.  From 
comparison  with  experiment  it  can  be  shown  that  the  model  underpredicts  the  ablation  mass  by 
about  25%.  However,  it  should  be  noted  that  some  mass  can  be  ablated  in  the  form  of  large 
particulates.  This  effect  for  one  particular  PPT  was  estimated  to  be  up  to  40%  of  the  total  ablated 
mass^^.  The  ablation  in  the  particulate  phase  was  not  considered  in  the  present  paper. 

The  near-cathode  sheath  model  allows  of  calculation  the  potential  drop  between  the  plasma  at 
the  thruster  exit  plane  and  the  cathode.  The  plasma  has  a  positive  potential  with  respect  to  the 
cathode  which  decreases  initially  with  time.  This  happens  because  the  current  pulse  is  essentially 
gone  after  5  its  when  the  main  plasma  cloud  arrives  at  the  exit  plane.  Therefore  the  plasma 
density  at  the  exit  plane  increases  while  the  current  decreases  and  this  leads  to  decreasing  A^. 

As  mentioned  earlier,  one  of  the  reasons  for  development  of  the  ID  cavity  model  Is  to  produce 
more  accurate  boundary  conditions  for  the  plume  simulation.  The  temporal  variation  of  the 
plasma  density  at  the  thruster  exit  plane  is  shown  in  Fig.  9.  In  comparison  with  previous  results, 
the  plasma  density  significantly  increases  at  later  times  since  plasma  parameter  variation  along 
the  cavity  length  is  taken  into  account.  These  new  boundary  conditions  are  employed  in  the 
particle  simulation  of  the  plume.  As  expected,  there  is  an  effect  on  the  plume  structure  at  later 
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times  as  shown  in  Fig.  10.  It  can  be  seen  that  the  simulation  predicts  well  the  initial  rise  and  the 
actual  peak  of  the  potential  data.  To  generate  the  results  shown  in  Fig.  10  we  assumed  that  the 
potential  at  the  thruster  exit  plane  is  constant  throughout  the  simulation.  However,  it  was  shown 
(Fig.  9)  that  this  potential  varies  with  respect  to  the  cathode.  The  measurements  of  plasma 
potential  were  made  with  respect  to  the  cathode^'^^  .  To  illustrate  this  effect  we  combine  the 
plasma  potential  in  the  plume  (Fig.  10)  with  variation  of  the  plasma  potential  at  the  thruster  exit 
plane  (Fig.  8).  These  results  are  shown  in  Fig.11  where  the  electron  temperature  near  the 
cathode  is  used  as  a  parameter.  One  can  see  that  this  approach  is  able  to  predict  the  drop  of  the 
potential  from  5  ps  and  the  minimum  before  the  rising  part.  In  making  this  comparison  with  the 
experimental  data,  it  is  expected  that  the  measurements  collected  over  the  first  5  ps  (before  mam 
plasma  cloud  arrives  at  the  thruster  exit)  cannot  be  reproduced  in  the  frame  of  the  present  model. 
This  is  due  to  dependence  on  the  ignition  plasma  introduced  by  the  igniter  spark,  which  is  not 
modeled  in  present  work.  In  the  present  model  we  assume  that  all  plasma  parameters  are  radially 
uniform  while  really  the  electron  temperature  near  the  nozzle  edges  may  be  smaller  than  that  at 
the  axis^®.  Therefore  we  introduce  the  electron  temperature  as  a  parameter. 


The  new  thruster  exit  boundary  conditions  provided  by  the  physical  models  presented  in  this 
study  do  not  significantly  affect  the  main  aspects  of  the  plume  structure  that  were  reported  in  Ref. 
7.  It  is  found  that  the  overall  distribution  of  the  chemical  species  in  the  plume  is  almost 
unchanged.  The  main  conclusion  from  Ref.  7  still  holds,  that  for  the  PPT-4,  the  main  component 
of  back  flow  onto  the  spacecraft  will  arise  from  the  highly  mobile  carbon  ions. 

Summary 

An  end-to-end  device/plume  model  has  been  developed  to  describe  plasma  generation, 
acceleration  and  plume  expansion  for  an  electrothermal  pulsed  plasma  thruster.  As  an  example, 
we  have  studied  the  Teflon  fed  PPT-4  developed  at  the  University  of  Illinois.  The  modeling 
strategy  was  based  on  a  one-dimensional  fluid  unsteady  model  for  the  plasma  generation  and 
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acceleration  processes  and  a  particle  simulation  of  the  plume.  A  new  kinetic  ablation  model  was 
employed  to  calculate  the  Teflon  ablation  rate  as  a  function  of  the  plasma  parameters.  Predicted 
results  were  compared  directly  with  data  for  PPT-4  performance  characteristics  and  plume 
flowfield.  In  general,  the  agreement  between  data  sets  was  good. 
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Figure  2.  The  ablation  rate  as  a  function  of  plasma  density  and  Teflon  surface  temperature.  The  arrow 
indicates  the  direction  of  ablation  rate  evolution  during  the  pulse. 
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Figure  3.  Variation  with  time  of  electron  temperature  and  ionization  degree  in  the  cavity 


stI  ‘auiix 


Figure  4.  Temporal  and  spatial  variation  of  the  chemical  species  during  the  discharge  pulse 
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Figure  7.  Ablated  mass  as  a  function  of  cavity  geometry  and  comparison  with  experiment,  a)  constant 
cavity  radius  Ra»  b)  constant  cavity  length  L. 
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Figure  8.  Temporal  variation  of  the  plasma  potential  at  the  thruster  exit  plane  with  respect  to  the  cathode 


Figure  9.  Variation  with  time  of  electron  density  and  comparison  with  previous  results 
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Figure  10.  Variation  of  the  plasma  potential  on  the  plume  centerline  at  10  cm  from  the  thruster  exit  plane 
and  comparison  with  experiment 
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Figure  11.  Variation  with  time  of  the  plasma  potential  in  the  plume  at  10  cm  from  the  PPT-4  exit  plane 
with  electron  temperature  near  the  cathode  as  a  parameter 
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The  vaporization  of  condensed  materials  in  contact  with  high-current  discharge  plasmas  is 
considered.  A  kinetic  numerical  method  named  direct  simulation  Monte  Carlo  (DSMC)  and 
analytical  kinetic  approaches  based  on  the  bimodal  distribution  function  approximation  are 
employed.  The  solution  of  the  kinetic  layer  problem  depends  upon  the  velocity  at  the  outer  boundary 
of  the  kinetic  layer  which  varies  from  very  small,  corresponding  to  the  high-density  plasma  near  the 
evaporated  surface,  up  to  the  sound  speed,  corresponding  to  evaporation  into  vacuum.  The  heavy 
particles  density  and  temperature  at  the  kinetic  and  hydrodynamic  layer  interface  were  obtained  by 
the  analytical  method  while  DSMC  calculation  makes  it  possible  to  obtain  the  evolution  of  the 
particle  distribution  function  within  the  kinetic  layer  and  the  layer  thickness.  ©  2001  American 
Institute  of  Physics.  [DOI:  10.1063/1.1345860] 


I.  INTRODUCTION 

The  vaporization  of  a  heated  surface  interacting  with  dis¬ 
charge  plasmas  has  a  great  interest  for  different  applications 
such  as  ablation  controlled  arcs,*  pulsed  plasma  thrusters, 
high-pressure  discharges,"^  vacuum  arcs,^  electroguns, ^  and 
metal  evaporation  by  laser  radiation  action.^  In  most  models, 
the  rate  of  evaporation  is  calculated  using  the  Langmuir 
relationship^  that  is,  however,  limited  to  the  case  of  vapor¬ 
ization  into  vacuum. 

Anisimov^  considered  a  case  of  vaporization  of  a  metal 
exposed  to  laser  radiation  using  a  bimodal  velocity  distribu¬ 
tion  function  in  the  nonequilibrium  (kinetic)  layer.  The  main 
result  of  this  work  is  the  calculation  of  the  maximal  flux  of 
returned  atoms  to  the  surface,  which  was  found  to  be  about 
18%  of  the  flux  of  vaporized  atoms.  This  result  was  obtained 
under  the  assumption  that  the  atom  flow  velocity  is  equal  to 
the  sound  velocity  at  the  external  boundary  of  the  kinetic 
layer.  In  many  physical  situations,  however,  the  expansion  of 
the  vapor  is  not  by  the  sound  speed  since  there  is  a  dense 
plasma  in  the  volume  discharge.  Beilis***’**  analyzed  metal 
vaporization  into  discharge  plasmas  in  the  case  of  a  vacuum 
arc  cathode  spot.  He  concluded  that  the  parameters  at  the 
outer  boundary  of  the  kinetic  layer  are  close  to  their  equilib¬ 
rium  values  and  that  the  velocity  at  the  outer  boundary  of  the 
kinetic  layer  is  much  smaller  than  the  sound  velocity.  In  both 
the  analyses  mentioned  above,  no  information  is  provided 
about  the  change  of  the  particle  velocity  distribution  function 
from  a  nonequilibrium  state  to  an  equilibrium  state  inside  the 
kinetic  layer. 

In  the  present  article  we  study  the  nonequilibrium  layer 
close  to  the  evaporating  surface  using  the  particle  method 
known  as  direct  simulation  Monte  Carlo  (DSMC).*^  It  will 
determine  the  thickness  of  the  nonequilibrium  layer  and  the 
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evolution  of  the  particle  distribution  function  within  the 
layer.  The  numerical  simulation  results  will  be  compared 
with  the  analytical  results  for  the  case  when  the  vapor  veloc¬ 
ity  at  the  outer  boundary  of  the  kinetic  layer  is  given  as  a 
parameter. 


II.  MODEL  OF  THE  NONEQUILIBRIUM  LAYER 

In  this  section  we  will  present  two  different  kinetic  mod¬ 
els  (particle  simulation  and  analytical  approach)  for  the  non¬ 
equilibrium  layer  near  the  evaporating  surface. 

A.  Particle  simulation 

In  the  nonequilibrium  layer  near  the  surface  there  are 
collisions  between  particles  that  eventually  lead  to  a  change 
of  the  distribution  function.  The  DSMC  method  uses  particle 
motion  and  collisions  to  perform  a  simulation  of  gas  dynam¬ 
ics  under  nonequilibrium  conditions.  Each  particle  has  spa¬ 
tial  and  velocity  coordinates.  The  collision  approach  between 
particles  is  based  on  a  probability  model  developed  from  the 
kinetic  theory  and  commonly  used  in  DSMC.*^ 

To  perform  the  DSMC  simulation  we  have  to  specify 
conditions  at  two  boundaries  (see  Fig.  1).  At  the  evaporating 
surface  with  density  Wq  temperature  Tq,  the  velocity 
distribution  function  for  emitted  particles  is  in  the  equilib¬ 
rium  form^’*** 


/o(V)  =  no 


l^rkTc] 


3/2 


expi 


mV^ 

2kfc 


V,>0. 


(1) 


At  the  outer  boundary  of  the  kinetic  layer  the  distribution 
function  for  particles  is  assumed  to  be 


/i(V)  =  ni 


^27TkT 


3/2 


expl 


(v,-v,y^vl+v: 


2kT^ 
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FIG.  1.  Schematic  representation  of  the  near  surface  layer. 


where  Vi  is  the  velocity,  is  the  density,  and  is  the 
temperature.  These  boundary  conditions  are  supplemented 
by  using  an  empirical  relation  between  Tq  and  /Zq.  As  an 
example,  we  have  used  the  equilibrium  vapor  pressure  for 
the  case  of  Teflon  in  the  form^ 

Pq  =  Pc^\^{-TcITq)  and  nQ-PQlkT^,  (3) 

where  Fq  is  the  equilibrium  pressure,  1.84X  lO^^N/m^ 
and  7^=20  815  K  are  the  characteristic  pressure  and  tem¬ 
perature  obtained  empirically.  The  physical  meaning  of  the 
characteristic  pressure  and  temperature  is  that  the  equilib¬ 
rium  pressure  equals  P^  when  vapor  reaches  the  temperature 
of  Tc .  The  calculations  are  performed  assuming  that  vapor 
consists  of  CF2  molecules  at  a  surface  temperature  in  the 
range  550-650  K  that  is  typical  for  an  electrothermal  pulsed 
plasma  thruster.^ 

The  DSMC  model  employed  in  the  present  article  has 
the  following  strategy.  Uniform  cells  0.5  X  in  size  are  em¬ 
ployed,  and  time  step  Ar  =  0.3X/y;„^j  (Ref.  13)  where  X  and 

_y=(2)tro/m)°'^  are  the  molecular  mean  free  path  and  the 
most  probable  thermal  speed  at  the  ablated  surface.  Mol¬ 
ecules  enter  the  flow  field  successively  from  the  surface  due 
to  evaporation,  and  from  the  outer  boundary  due  to  Maxwell¬ 
ian  velocity  distribution  that  allows  particle  velocities  in  the 
negative  direction.  Using  the  assumption  about  the  distribu¬ 
tion  function  at  the  surface  and  at  the  external  boundary  Eqs. 
(1)  and  (2)  we  can  calculate  the  flux  of  molecules  entering 
from  the  surface,  ,  and  from  the  outer  boundary  of  the 
layer,  G^ ,  as  follows: 

G,=no(2kTolm)°-^l2'iT'^-^  (4) 

and 

Gi,  =  n^{2kTi 

X  [exp(  -  a^)  -  1  -  erf(  a)}]/2'7r®-^,  (5) 

where  a-V^  l{2kTi  The  molecular  interaction  is  de¬ 

scribed  by  the  variable  hard-sphere  (VHS)  model. The 
VHS  model  is  employed  to  select  molecular  collision  pairs 
from  cells  and  to  distribute  the  postcollision  velocities.  This 
model  assumes  that  the  scattering  from  molecular  collision  is 
isotropic  in  the  center  of  mass  frame  of  reference. Both 
boundaries  (wall  and  outer  boundary  of  the  kinetic  layer)  are 
assumed  to  be  perfectly  absorbing.  The  flow  will  arrive  at  a 


steady  state  when  the  sum  of  G^  and  G^  is  exactly  balanced 
by  the  flux  of  molecules  leaving  from  the  outer  boundary  or 
sticking  on  the  surface: 


G,+  G£,  =  G,+G^,  (6) 

where  G^  is  the  flux  of  the  particles  returned  to  the  surface 
during  the  time  step  and  Gy  is  the  flux  of  the  particles  cross¬ 
ing  the  outer  boundary  of  the  layer. 

In  the  DSMC  approach,  in  order  to  calculate  the  evolu¬ 
tion  of  the  distribution  function  inside  the  kinetic  layer,  we 
have  to  specify  the  thickness  of  the  layer  in  units  of  X.  The 
parameters  at  the  outer  boundary  of  the  kinetic  layer  (ni  and 
Tx)  and  the  flux  of  returned  particles  are  calculated  as  a 
function  of  the  distance  of  the  location  of  the  outer  boundary 
of  the  layer  and  of  the  velocity  at  this  boundary  Vj . 


B.  Analytical  approach 

Let  us  consider  the  analytical  approach  developed  in 
Refs.  9-11,  where  the  vapor  parameters  Tj  and  nj  at  the 
outer  boundary  can  be  obtained  without  information  about 
the  layer  thickness.  This  means  that  the  problem  is  reduced 
to  integration  of  the  conservation  equations  of  mass,  momen¬ 
tum,  and  energy  across  the  layer.^  We  consider  a  nonequi¬ 
librium  layer  (thickness  of  about  a  mean  free  path  X)  adja¬ 
cent  to  the  surface  (as  shown  in  Fig.  1),  where  the  velocity 
distribution  function  of  the  evaporated  molecules  reached 
equilibrium  by  the  rare-field  collisions  with  the  background 
heavy  particles  and  furthermore  the  vapor  flow  is  described 
by  a  hydrodynamic  approach.  Using  Anisimov’s  assumption^ 
that  the  velocity  distribution  function  for  the  returned  par¬ 
ticles  ( Vx^O)  is  /3/i(V),  where  p  is  the  proportionality  co¬ 
efficient,  the  relation  of  the  heavy  particle  parameters  at  the 
outer  boundary  of  the  kinetic  layer  in  the  case  of  an  arbitrary 
velocity  is  obtained  from  the  model and  reads  as  follows: 

no  ^ 

X  {exp(  —  a^)  —  erfc(  a)}. 


no 

4^ 


-  /3[(0.5+  a^)erfc(a)  -  a  exp(  - 


(7) 


X  erfc(  a)  -  (2  +  a^)exp  ( -  a'^)/'7r®’^}], 

where  dQ  —  mf2kTQ,  di=mf2kTx ,  erfc(a)=l— erf(a),  and 
erf(a)  is  the  error  function.  The  equation  system  (7)  is  ob¬ 
tained  using  the  boundary  conditions  (l)-(3)  and  the  conser¬ 
vation  laws  of  mass,  momentum,  and  energy  across  the 
layer.^’^®  By  calculating  the  parameters  at  the  outer  boundary 
of  the  kinetic  layer  we  can  obtain  the  flux  of  returned 
particles: 
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FIG.  2.  Comparison  of  the  analytic  and  DSMC  return  flux  as  a  function  of 
velocity  Vj  with  the  distance  of  location  of  the  outer  boundary  of  the  kinetic 
layer  L  as  a  parameter. 


FIG.  3.  The  DSMC  calculated  return  flux  as  a  function  of  the  distance  of 
location  of  the  outer  boundary  of  the  kinetic  layer  L  in  the  case  of  sound 
velocity  at  the  boundary  1. 


/: 


pf,{V)VdV 


I  kT 

=  T — ^  {exp(-a^)-a7r®-^erfc(a)}.  (8) 

\Z7Tm  I 

The  system  of  equations  (7)  has  four  unknowns  and  therefore 
the  solution  can  be  found  having  one  unknown  as  a  param¬ 
eter,  which  is  the  velocity  Vj  at  the  outer  boundary  of  the 
kinetic  layer  in  our  case. 


III.  RESULTS 

DSMC  calculations  and  a  comparison  with  the  analytical 
predictions  for  the  flux  of  returned  atoms  is  shown  in  Fig.  2 
where  the  thickness  of  the  kinetic  layer  is  used  as  a  param¬ 
eter.  One  can  see  that  in  the  case  of  small  velocity  (a 
0.5)  at  the  outer  boundary  all  results  agree  well.  This  is  the 
case  when  the  thickness  of  the  nonequilibrium  layer  is  about 
one  mean  free  path.  In  the  case  when  evaporation  occurs  at 
about  the  sound  velocity  at  the  outer  boundary,  the  DSMC 
calculations  approach  the  analytical  value  at  a  layer  thickness 
of  '-'10-20  mean  free  paths. 

The  calculation  of  the  backflux  dependence  with  dis¬ 
tance  inside  the  layer  when  a  =  1  is  presented  in  Fig.  3.  It 
can  be  seen  that  in  the  case  when  Vj  is  about  the  sound 
velocity,  the  flux  of  the  returned  molecules  depends  upon  the 
distance  from  the  evaporating  surface  where  the  external 
boundary  is  placed.  Thus  up  to  a  layer  thickness  of  about  20 
mean  free  paths,  the  flux  changed  strongly  and  then  it  is 
weakly  saturated.  The  DSMC  calculation  predicts  a  16%  flux 
of  returned  particles,  which  is  very  close  to  the  analytical 
result  of  18%.  The  reason  for  this  difference  can  be  under¬ 
stood  by  analyzing  the  velocity  distribution  function  of  re¬ 
turned  particles  in  the  DSMC  calculation. 

Results  of  the  DSMC  calculation  of  the  velocity  distri¬ 
bution  function  and  comparison  with  the  analytic  approxima¬ 
tion  are  shown  in  Fig.  4(a)  for  the  case  of  sound 

velocity  at  the  outer  boundary.  One  can  see  that  the  distribu¬ 
tion  functions  remain  different  in  the  case  of  a  lOOX  kinetic 


layer  thickness.  This  is  not  the  case  when  the  velocity  at  the 
outer  boundary  of  the  kinetic  layer  is  small  as  shown  in  Fig. 
4(b),  where  the  DSMC  distribution  function  agrees  well  with 
the  analytic  approximation.  Therefore  it  is  not  surprising  that 
the  calculated  flux  of  returned  particles  is  also  found  to  be  in 
good  agreement  with  the  analytical  result.  It  should  be  noted 


(a)  V/(2kT»°" 


(b)  V/(2kT^m)°* 


FIG.  4.  Variation  of  the  velocity  distribution  function  of  the  returned  par¬ 
ticles  near  the  wall  with  the  distance  of  location  of  the  outer  boundary  of 
the  kinetic  layer  L  as  a  parameter,  (a)  Vi  =  {5kTi/3m)^^  and  (b)  Vj 
=  0.1(2itri  /mf\ 
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FIG.  5.  Variation  of  the  velocity  distribution  function  (normal  to  the  surface 
component)  with  the  distance  from  the  evaporated  surface  as  a  parameter. 
The  thickness  of  the  kinetic  layer  lOOX  and  Vi  =  (5itrj/3m)°  ^ 

that  the  discontinuity  in  the  analytical  distribution  function 
[Fig.  4(b)]  is  the  result  of  the  definition  of  the  distribution 
function  of  the  returned  particles  (Sec.  II B). 

The  evolution  of  the  particle  distribution  function  within 
the  Knudsen  layer  is  shown  in  Fig.  5  for  the  case 
=(5kTil3m)^'^.  One  can  see  that  the  velocity  distribution 
function  approaches  a  drifted  Maxwellian  at  a  distance  of 
several  mean  free  paths  from  the  surface.  The  drift  velocity 
slightly  increases  with  further  distance  from  the  evaporating 
surface. 

The  results  of  calculation  of  the  analytic  system  of  equa¬ 
tions  (7)  are  presented  in  Fig.  6  with  the  normalized  velocity 
Vj  as  a  parameter.  The  temperature  Tj ,  density  nj ,  and  the 
flux  of  returned  particles  J-  all  decrease  as  the  velocity  at 
the  outer  boundary  of  the  kinetic  layer  increases.  In  the  lim¬ 
iting  case  of  the  sound  velocity,  the  flux  of  returned  particles 
is  equal  to  18%  as  was  obtained  by  Anisimov.^  In  this  case 
the  analytically  predicted  flux  of  returned  particles  is  larger 
than  that  obtained  by  numerical  simulations  (16%,  see  Fig. 
3).  It  should  be  pointed  out  that  the  dependence  of  the  flux  of 


returned  particles  7  _  on  the  velocity  Vj  has  a  minimum  near 
the  sound  speed  (see  Fig.  6).  The  minimum  corresponds  to 
the  sound  speed  with  adiabatic  index  1.3.  This  fact  as  well  as 
the  overestimate  of  the  returned  particle  flux  is  connected 
with  the  assumption  in  the  analytical  approach  of  the  form  of 
the  particle  velocity  distribution  for  the  particles  returned  to 
the  evaporating  surface,  i.e.,  that  the  distribution  function  of 
the  returned  particles  is  proportional  to  the  distribution  func¬ 
tion  at  the  outer  boundary  of  the  kinetic  layer.^ 

IV.  SUMMARY 

Two  kinetic  approaches,  namely  the  particle  method 
(DSMC)  and  bimodal  distribution  function  approach  were 
employed  to  describe  the  parameters  in  the  nonequilibrium 
kinetic  layer  near  the  evaporating  surface.  DSMC  calculation 
makes  it  possible  to  establish  the  thickness  of  the  kinetic 
layer  and  the  evolution  of  the  particle  distribution  function 
within  the  layer.  The  thickness  of  the  kinetic  layer  and  the 
vapor  density  and  temperature  in  the  kinetic  layer  adjacent  to 
the  evaporating  surface  depend  upon  the  velocity  at  the  outer 
boundary  of  the  layer.  We  have  found  that  the  thickness  of 
the  kinetic  layer  increases  from  a  few  mean  free  paths  X  in 
the  case  of  small  velocity  up  to  about  10-20  X  in  the  case  of 
the  evaporation  with  sound  speed  at  the  outer  boundary  of 
the  layer.  Comparison  of  the  DSMC  and  analytical  results 
indicates  that  the  analytical  model  predicts  correctly  the  flux 
of  returned  particles  over  a  wide  range  of  velocity  at  the 
outer  boundary  of  the  layer.  The  present  model  can  be  used 
for  calculation  of  the  rate  of  evaporation  of  the  heated  sur¬ 
face  interacting  with  a  plasma.  The  free  parameter  of  this 
model,  the  velocity  at  the  outer  boundary  of  the  layer,  can  be 
determined  by  coupling  this  model  with  a  model  of  the  hy¬ 
drodynamic  layer  and  the  plasma  bulk. 
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Abstract 

A  kinetic  model  is  developed  of  Teflon  ablation  caused  by  a  plasma. 

The  model  takes  into  account  the  returned  atom  flux  that  forms  in  the 
non-equilibrium  layer  during  the  ablation.  This  approach  makes  it  possible 
to  calculate  the  ablation  rate  for  the  case  when  the  Teflon  surface 
temperature  and  the  density  and  temperature  in  the  plasma  bulk  are  known. 


The  problem  of  the  ablation  controlled  discharge  has  a  common 
general  interest  in  a  number  of  applications  such  as  electric 
fuses,  circuit  breakers,  soft  x-ray,  pulsed  plasma  thrusters 
and  extreme  ultraviolet  sources  [1-5].  In  these  devices, 
the  discharge  energy  is  principally  dissipated  by  ablation  of 
wall  material,  which  then  forms  the  main  component  of  the 
discharge  plasma.  The  ablated  vapour  increases  the  pressure 
within  the  capillary  and  the  plasma  is  expelled  through  the  exit. 

Previously,  most  of  the  plasma  models  of  the  ablated 
controlled  discharge  employed  the  Langmuir  relationship  [6], 
which  is  limited  to  the  case  of  material  ablation  into  a  vacuum. 
For  the  conditions  of  a  pulsed  plasma  thruster  (PPT)  this 
approach  was  also  recently  used  [7, 8],  However,  the  process 
of  ablation  in  the  ablation  controlled  discharge  should  be 
described  taking  into  account  the  fact  that  in  the  Teflon  cavity 
the  vapour  does  not  expand  into  vacuum  but  rather  into  a 
volume  discharge.  For  application  to  the  vacuum  arc  discharge, 
a  kinetic  model  of  the  cathode  vapourization  in  the  non¬ 
equilibrium  plasma  layer  was  developed  by  Beilis  [9, 10].  It 
was  shown  that  the  flux  of  returned  atoms  toward  the  surface 
can  become  comparable  to  the  flux  of  vaporizing  atoms.  It  was 
also  concluded  [9, 10]  that  the  parameters  at  the  outer  boundary 
of  the  kinetic  layer  are  close  to  their  equilibrium  values  and  that 
the  velocity  at  the  outer  boundary  of  the  kinetic  layer  is  much 
smaller  than  the  sound  velocity.  Therefore,  it  was  found  that 
the  cathode  erosion  rate  in  an  arc  discharge  by  vapourization  is 
much  smaller  than  the  solid  body  evaporation  rate  into  vacuum. 
In  the  present  work  we  employ  a  kinetic  approach  similar  to 
that  used  for  the  cathode  vacuum  arc  evaporation  to  calculate 
Teflon  ablation  parameters.  As  an  example,  the  conditions 
typical  for  an  electrothermal  PPT  are  considered  [4, 7]. 

The  problem  starts  by  considering  the  non-equilibrium 
kinetic  layer  near  the  evaporating  surface.  Let  us  consider 
the  structure  of  the  near  surface  layers  in  detail  as  shown  in 
figure  1.  One  can  distinguish  two  different  layers  between 
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Figure  1,  Schematic  representation  of  the  layer  structure  near  the 
ablated  surface. 


the  surface  and  the  plasma  bulk:  (1)  a  kinetic  non-equilibrium 
layer  adjusted  to  the  surface  with  a  thickness  of  a  few  mean 
free  paths  (the  Knudsen  layer),  (2)  a  collision-dominated 
non-equilibrium  layer,  where  the  electron  and  heavy  particle 
temperatures  differ.  It  is  assumed  that  at  the  right  edge  of  the 
second  layer,  all  species  (ions,  neutrals  and  electrons)  reach 
thermal  equilibrium.  The  basic  idea  of  the  present  model 
consists  in  combining  the  model  for  the  kinetic  layer  [9-11] 
with  hydrodynamic  layer  analyses  in  the  second  layer.  Use 
of  the  general  plasma  model  [7]  in  the  equilibrium  region 
provides  the  electron  temperature  T2  and  plasma  density  «2- 
Firstly,  we  briefly  introduce  the  kinetic  model  of  the  non¬ 
equilibrium  layer.  According  to  the  approach  of  the  work 
in  [9]  and  [10]  (using  Anisimov’s  [11]  assumption  that  the 
velocity  distribution  function  for  the  returned  particles  is 
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where  is  a  proportionality  coefficient  and  fi(V) 
is  the  Maxwellian  distribution  function  shifted  by  Vi,V  is  the 
velocity  vector)  the  relation  of  the  heavy  particle  parameters 
at  the  outer  boundary  of  the  kinetic  layer  in  the  case  of 
arbitrary  velocity  Vi  are  as  described  by  the  following  set  of 
equations: 


^  =  :^((l  +2«^)-;6[<0-5+a^)erfc(a) 

4do  2d\ 

—a  exp(— a^)/7r°'^) 


no 


11  =  ^^■'[“(“^+2.5) 

(^/2){(2.5+aVerfc(a)  -  (2+a2)exp(-a^)/;r°-^] 


(1) 


where  4)  =  mllkTo,  d\  =  ml2kT\,  a  =  V\l{2kT\lrnf-^, 
erfc(a)  =  1  erf  (a),  erf  (or)  is  the  error  function, 
Tq  is  the  surface  temperature  and  no  is  the  equilibrium 
density. 

Very  recently  it  was  shown  that  the  velocity  at  the  outer 
boundary  of  the  kinetic  layer  Vi  strongly  affects  the  kinetic 
layer  parameters  [12].  To  find  the  velocity  Vi  we  apply 
the  mass  and  momentum  conservation  equations  for  heavy 
particles  in  the  hydrodynamic  region  (assuming  a  single  fluid 
model)  between  boundaries  2  and  3.  Assuming  weakly  ionized 
plasma  in  the  hydrodynamic  layer,  the  integration  of  the  mass 
and  momentum  conservation  equations  yields  the  following 
relations  between  parameters  at  boundaries  2  and  3: 


mV,  =n2V2  (2) 

nikTi+mniVf  —  n2kT2  +  mn2V2‘  (3) 

Combination  of  these  two  equations  yields  the  following 
expression  for  the  velocity  at  the  outer  boundary  of  the  kinetic 
(Knudsen)  layer: 


V^/(2kTi/m)  =  (T2n2/2Ti  -  «i/2)/(n,  -  n]/n2).  (4) 


This  equation  makes  it  possible  to  calculate  the  velocity  at 
boundary  1  and  therefore  to  calculate  the  ablation  rate  that  is 
proportional  to  Vin, .  The  system  of  equations  is  closed  if  the 
equilibrium  vapour  pressure  can  be  specified.  In  the  case  of 
Teflon,  the  equilibrium  pressure  formula  is  used: 

P  =  P,exp(-Tc/To)  (5) 

where  P  —  nokTo  is  the  equilibrium  pressure.  Pc  —  1.84  x 
10^^  N  m“^  and  7^  =  20  815  K  are  the  characteristic  pressure 
and  temperature,  respectively  [4]. 

The  solution  of  the  problem  depends  upon  plasma  density 
^2,  plasma  temperature  T2  and  surface  temperature  Tq.  The 
parameters  «2»  T2  are  determined  by  the  bulk  plasma  flow.  It 
was  estimated  from  various  experiments  that,  under  typical 
PPT  operation  conditions,  the  plasma  density  near  the  Teflon 
surface  is  about  lO^’-lO^"^  m"^,  the  plasma  temperature  is 
about  1-4  eV  and  the  Teflon  surface  temperature  Tq  is  about 
600-650  K  [4,7,8,13].  In  the  present  paper  we  present 
solutions  with  n2,  T2  and  Tq  as  parameters  in  the  ranges 
mentioned  above. 


Figure  2,  The  velocity  V,  as  a  function  of  plasma  temperature  with 
Teflon  surface  temperature  as  a  parameter. 


Figure  3.  The  ablation  rate  as  a  function  of  plasma  temperature 
with  Teflon  surface  temperature  as  a  parameter. 

The  dependence  of  the  velocity  V,  on  the  electron 
temperature  T2  is  shown  in  figure  2  with  surface  temperature 
Tq  as  a  parameter  for  given  7x2 .  One  can  see  that  the  velocity  V, 
remains  small  over  the  entire  range  of  plasma  temperature  and 
generally  decreases  with  temperature  increase.  The  velocity 
V,  is  very  sensitive  to  the  plasma  temperature  variation  in  the 
case  of  relatively  small  surface  temperature.  As  a  result  of  this 
dependence,  the  ablation  rate  also  decreases  with  increasing 
electron  temperature  as  shown  in  figure  3. 

Ablation  rate  contours  in  the  plane  with  the  plasma  density 
and  Teflon  surface  temperature  as  the  coordinates  are  displayed 
in  figure  4.  The  same  ablation  rate  that  can  be  found  in  the 
high  and  low  density  range  corresponds  to  the  solution  of  the 
problem  with  small  and  large  velocity  at  the  outer  boundary 
of  the  kinetic  layer,  respectively.  There  is  no  solution  for  the 
ablation  rate  in  regions  above  the  curve  with  ablation  rate  equal 
to  zero.  This  region  in  the  712— To  plane  corresponds  to  the 
case  when  the  right-hand  side  in  equation  (4)  is  negative.  The 
physical  meaning  of  these  results  can  be  explained  as  follows. 
In  the  ablation-controlled  discharge,  the  plasma  density  in  the 
bulk  is  determined  by  the  rate  of  ablation  from  the  surface. 
In  the  case  of  small  surface  temperature  one  can  expect  a 
smaller  ablation  rate  and  therefore  high  plasma  densities  in 
the  discharge  cannot  be  generated. 
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Teflon  surface  temperature.  K 

Figure  4.  Contours  of  ablation  rate  in  the  plasma  density 
(n2)-Teflon  surface  temperature  (Tq)  plane. 


24  kg  s“^  which  is  close  to  that  measured  in  the 
experiment  (~29  kg  m”^  s~^  [13]). 

In  summary,  we  have  developed  a  kinetic  model  of  the 
material  ablation  in  an  ablated  controlled  discharge  with 
application  to  a  pulsed  plasma  thruster.  The  model  accounts  for 
the  case  when  the  velocity  at  the  outer  boundary  of  the  kinetic 
layer  is  smaller  than  the  sound  velocity  due  to  the  presence 
of  a  high-density  discharge  plasma.  The  present  model  can 
be  coupled  with  a  plasma  discharge  model  to  describe  the 
electrical  discharge  self-consistently. 
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[4, 13].  During  the  discharge  pulse,  the  plasma  parameters 
vary  so  that  near  the  current  peak  one  can  expect  an  ablation 
rate  higher  than  average  while  towards  the  end  of  the  pulse 
the  ablation  rate  decreases.  The  time-averaged  ablated  mass 
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Model  of  Particulate  Interaction  with  Plasma 
in  a  Teflon  Pulsed  Plasma  Thruster 

Michael  Keidar*  and  Iain  D.  Boyd'*^ 

University  of  Michigan,  Ann  Arbor,  Michigan  48109 
and 

Isak  L  Beilis^ 

Tel  Aviv  University,  Tel  Aviv  69978,  Israel 

The  presence  of  particulates  [referred  to  as  macroparticles  (MPs)]  in  the  plume  of  a  pulsed  plasma  thruster  (PPT) 
may  affect  interaction  with  the  spacecraft.  Possible  particulate  related  effects  depend  on  particulate  properties. 
The  MPs  emitted  into  the  plasma  during  the  discharge  may  be  charged,  accelerated,  and  heated  by  the  ion, 
electron,  and  neutral  fluxes  depending  on  the  MP  residual  time.  Different  aspects  of  MP-plasma  interaction  in 
the  experimentally  observed  range  of  MP  radii  (0.1-100  ^m)  are  analyzed.  It  is  found  that  the  chaining  time  is 
smaller  whereas  the  steady-state  potential  is  larger  in  the  case  of  a  large  MP.  A  1-^m  MP  is  found  to  have  a  charge 
of  about  10^  electrons  in  the  case  of  an  electron  density  of  about  10^  m“^.  The  two  primary  forces  acting  on  the 
MP  in  the  PPT  discharge  are  a  drag  force  due  to  collision  with  neutral  atoms  and  ions  and  an  electric  force  due 
to  the  presence  of  the  electric  field  in  the  current  cariying  plasma.  The  calculation  of  the  MP  velocity  shows  that 
the  maximum  possible  velodty  of  a  l-yitm  MP  is  about  230  m/s,  which  is  close  to  that  estimated  experimentally. 
Only  small  MPs  (~  0.1  /xm)  can  be  entrained  by  the  plasma  jet,  whereas  large  MPs  are  generally  slower  and  flow 
substantially  behind  the  plasma  jet.  MP  temperature  and  decomposition  rates  are  calculated  by  solving  a  heat 
balance  equation.  It  is  found  that  small  MPs  (<1  ^m)  may  completely  decompose  during  a  1-^s  pulse. 


Nomenclature 

C  =  specific  heat 

Cs  =  sound  speed 

E  =  electric  field 

Fd  =  total  drag  force 

g  =  Teflon®  ablation  rate 

/  =  discharge  current 

Je{r)  =  electron  current  density  in  the  sheath  around 
microparticle  (MP) 

JeO  =  electron  thermal  flux  at  the  plasma- sheath  interface 
Ji{r)  =  ion  current  density  in  the  sheath  around  MP 
7,0  =  ion  flux  in  the  absence  of  a  field  at  the  plasma-sheath 

interface 

j  =  current  density 

L  =  cavity  length 

Ld  -  Debye  length 

Mp  =  MP  mass 

=  atom  mass 

N  =  plasma  density  at  the  plasma- sheath  interface 

Na  =  plasma  density  near  the  anode,  that  is,  where  z  =  0 

Nu  =  Nusselt  number 

n  =  plasma  density  normalized  by  Na 

Pq  =  equilibrium  pressure 

Qd  =  cooling  rate  due  to  decomposition  of  material 

Qi  e  =  heat  rate  due  to  ion  and  electron  flux 

Q„  =  neutral  atom  heat  flux 

Qp  =  MP  charge 

Qr  =  radiation  cooling  rate 

Rp  =  MP  radius 

Rpo  =  initial  MP  radius 
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Tfl  =  radius  of  cavity 

T  =  neutral  gas  temperature 

T^  -  electron  temperature 

Ti  =  ion  temperature 

Tp  =  MP  temperature 

Up  -  MP  potential  with  respect  to  the  plasma 

Us  =  potential  drop  in  the  near  spacecraft  surface  sheath 

V  =  plasma  velocity 

Vd  =  MP  ablation  rate 

Vp  =  MP  velocity 

Vp„  =  normal  component  of  the  MP  velocity 
u  =  plasma  velocity  normalized  by  the  sound  speed 
z  =  coordinate  in  the  axial  direction  normalized  by  cavity 
length 

a  =  heat  transfer  coefficient 

AH  =  ablation  heat 

AT  =  temperature  difference  between  the  plasma  and  the  MP 

At  =  residual  time  of  MP  in  discharge 

s  =  dielectric  permittivity  of  Teflon® 

fo  =  permittivity  of  vacuum 

kc  =  mean  free  path 

I  =  particle  emissivity 

p  =  specific  density 

a  =  plasma  conductivity 

Or  =  Stefan-Boltzmann  constant 

r  =  dimensionless  time,  Ld/C, 

(f)  =  normalized  MP  potential, 

I.  Introduction 

HE  pulsed  plasma  thruster  (PPT)  has  been  reconsidered  re¬ 
cently  as  an  attractive  spacecraft  propulsion  option.  This  has 

happened  mainly  due  to  a  greater  emphasis  being  placed  on  the 
development  of  satellites  with  reduced  size  for  many  applications. 
PPTs  are  expected  to  provide  exact  impulse  bits  to  be  used  for  accu¬ 
rate  attitude  control.  The  principal  advantage  of  PPTs  is  their  simple 
design,  which  provides  high  reliability.  In  particular,  the  higher  re¬ 
liability  of  the  PPT  is  achieved  through  the  use  of  solid  propellant, 
which  eliminates  design  and  operation  complexity  connected  with 
using  liquid  and  gas  propellants. 
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Unfortunately,  the  PPT  has  a  very  poor  performance  characteris¬ 
tic.  For  instance,  a  flight-qualified  PPT  design  had  an  efficiency  of 
about  8%  (Ref.  3).  One  of  the  factors  leading  to  low  efficiency  is  the 
late  neutral  ablation.'^  Another  factor  that  may  significantly  reduce 
the  efficiency  is  the  particulate  emission.  Estimates  have  shown  that 
the  particulate  emission  consumes  about  40%  of  the  total  propellant 
mass,  albeit  contributing  only  1%  to  the  total  thrust.^ 

Particulates,  sometimes  referred  to  as  macroparticles  (MPs),  are 
emitted  during  the  pulse  and  may  interact  with  the  surrounding 
plasma.  Various  experimental  and  theoretical  aspects  of  MP  inter¬ 
action  with  plasma  have  been  studied  in  diiferent  systems  such  as 
plasma  assistedchemical  vapor  deposition,^  plasma  etching^  rf  and 
glow  discharge cathodic  vacuum  arc  deposition,  ^  ^  and  in- space 
plasma.'^  The  MPs  are  subjected  to  charging,  heating,  and  momen¬ 
tum  transfer.^"’^  It  may  be  concluded  that  in  spite  of  many  common 
features  of  MP- plasma  interaction,  the  MP  charge,  velocity,  and 
temperature  depend  on  properties  of  the  plasma  and  MP  related  to 
each  specific  system.  Recently,  results  for  the  MP  size  distribution  in 
a  PPT  have  been  presented  by  Spanjers  et  al.^  Particulates  were  ob¬ 
served  to  have  characteristic  diameters  ranging  from  about  0.1  fxm 
to  over  1(X)  fim.  The  PPT  plasma  has  a  weak  degree  of  ionization.^ 
Thus,  one  can  expect  an  important  role  of  the  neutral  component 
in  the  momentum  and  energy  transfer  to  MPs.  No  analyses  of  MP 
interaction  with  the  PPT  plasma  have  been  reported  previously.  The 
present  work  has  the  primary  objective  of  removing  this  deficiency 
by  analyses  of  different  aspects  of  MP  behavior  of  a  PPT,  such  as 
charging,  acceleration,  heating,  and  decomposition.  This  work  will 
provide  additional  information  about  MPs  to  help  mission  planners 
in  estimating  the  particulate  contamination  potential. 

IL  Model  of  MP-Plasma  Interaction 

In  the  present  section  we  will  develop  the  model  of  MP  interaction 
with  plasma,  including  MP  charging,  transport,  and  heating  in  a 
pulsed  discharge.  As  a  working  example,  we  will  concentrate  on 
a  specific  type  of  PPT,  developed  at  the  University  of  Illinois,  the 
so-called  PPT-4  (Ref.  1).  This  PPT  has  a  coaxial  configuration  in 
which  Teflon®  is  ablated  from  a  cylindrical  cavity  sitting  in  front 
of  the  central  electrode  of  6-mm  diam  and  an  annular  electrode  of 
43 -mm  diam.  The  two  electrodes  are  connected  with  a  30-deg  half¬ 
angle  nozzle.  The  typical  PPT-4  pulse  duration  is  about  10  ixs  with 
a  current  peak  of  about  8  kA.  The  main  physical  process  in  this 
type  of  thruster  occurs  in  the  Teflon  cavity.  Rapid  heating  of  a  thin 
dielectric  surface  layer  leads  to  decomposition  of  the  material  of  the 
wall.  As  a  result  of  heating,  decomposition,  and  partial  ionization  of 
the  decomposition  products,  the  total  number  of  particles  increases 
in  the  cavity. 

At  the  same  time,  nonuniformities  in  the  discharge  distribution 
across  the  Teflon  surface  may  cause  overheating  followed  by  phase 
change  of  propellant.^  High  plasma  pressure  in  the  PPT  channel  may 
lead  to  Teflon  MP  generation.  However,  the  mechanism  of  such 
MP  generation  from  the  propellant  is  not  understood  sufficiently. 
Another  possible  source  of  MPs  is  the  spot  attachment  at  electrodes 
that  is  a  typical  phenomenon  in  the  early  stages  of  discharge.^^ 
Several  scenarios  leading  to  particle  generation  from  the  electrodes 
are  possible:  action  of  the  plasma  pressure  on  the  liquid  pool  in  the 
quasi-steady  regime  may  form  droplets,  similar  to  that  occurring  in 
the  vacuum  arc  cathode  spot^'*;  growth  of  the  perturbation  of  the 
liquid  surface  due  to  the  Rayleigh- Taylor  instability  may  result  in 
liquid  jet  fragmentation  and  small  droplet  generation,  as  it  occurs  in 
liquid- metal  ion  sources solid-particle  generation  due  to  a  wave 
of  thermoelastic  stresses  may  also  be  caused  by  local  overheating.^^ 

In  the  experiment,  basically  two  populations  of  MPs  were  found.  ^ 
One  population  of  MPs  is  characterizedby  a  diameter  less  than  1  ^m 
and  by  a  spherical  shape.  It  was  concluded  that  these  MPs  come  from 
the  steel  PPT  electrodes^  due  to  discharge  attachment.  The  second 
population  has  a  more  random  shape  and  size  with  maximal  size  up 
to  200  fim.  The  second  population  is  concluded  to  be  originated  from 
the  Teflon  propellant.^  In  PPT-4,  where  copper  electrodes  were  used, 
tracks  of  discharge  attachment  at  the  electrodes  were  not  observed, 
as  was  noted  in  a  private  communication  with  R.  L.  Burton.  How¬ 
ever,  Teflon  macroparticles  may  still  be  generated  that  have  a  sig¬ 
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Fig.  1  Schematic  presentation  of  MP-plasma  interaction. 


nificant  effect  on  propellant  losses.  Thus,  as  we  consider  the  PPT-4 
configuration  in  the  present  paper,  we  will  concentrate  on  these  MPs. 

The  MPs  emitted  into  the  plasma  during  the  discharge  may  be 
charged,  accelerated,  and  heated  by  the  ion,  electron,  and  neu¬ 
tral  fluxes,  depending  on  the  MP  residual  time.  It  is  assumed  in 
the  model  that  the  MP  may  be  emitted  at  any  point  of  time  during  to 
the  discharge  pulse.  This  implies  that  MPs  generatedat  the  laterstage 
of  discharge  will  have  less  time  for  interaction  with  plasma.  The 
scheme  of  MP-plasma  interaction  is  shown  in  Fig.  1.  The  present 
model  is  based  on  the  assumption  that  during  the  MP  flight  the 
plasma  parameters  do  not  change  significantly,  that  is,  the  follow¬ 
ing  condition  is  fulfilled:  VpAt  where  X  is  any 

of  the  plasma  parameters  (velocity,  density,  and  temperature).  This 
model  of  plasma- MP  interaction  is  approximation  to  the  more  com¬ 
plex  reality  in  which  plasma  parameters  vary  during  the  MP  residual 
time.  In  the  present  model,  we  will  analyze  the  plasma-MP  interac¬ 
tion  having  plasma  density  and  velocity  that  lie  in  the  range  of  their 
possible  variation  during  the  discharge  pulse,  as  free  parameters  of 
the  problem. 


A.  Quasi-Steady-State  Plasma  Model 

The  main  features  of  the  electrical  discharge  in  the  dielectric  cav¬ 
ity  include  joule  heating  of  the  plasma,  heat  transfer  to  the  Teflon,  de¬ 
composition  followed  by  ionization,  and  acceleration  of  the  plasma 
up  to  the  sound  speed  at  the  cavity  exit.  In  this  section,  we  will 
present  a  simple  quasi-stationary  model  of  the  plasma  flow  in  the 
Teflon  cavity.  A  quasi-steady-state  approach  to  the  PPT  flow  has 
numerous  limitations.’^  It  requires  that  the  propellant  ablation  must 
supply  material  to  the  discharge  chamber  in  times  shorter  than  the 
characteristic  flow  convection  time,  which,  in  turn,  should  be  less 
than  the  characteristic  time  of  discharge  parameter  variation.  How¬ 
ever,  this  approach  may  provide  some  useful  information  about  the 
possible  range  of  density  in  the  discharge  chamber  and  the  spatial 
plasma  density  and  velocity  distribution  along  the  cavity  length.  In 
the  present  model  the  plasma  velocity  and  density  will  be  used  in 
the  MP-plasma  interaction  model  as  parameters. 

We  apply  a  one-dimensional  hydrodynamic  model  for  the  plasma 
under  the  assumption  that  the  Teflon  evaporates  uniformly.  Products 
of  Teflon  evaporation  are  partially  ionized  in  the  cavity.  Partially 
ionized  plasma  conducts  the  current,  which  is  carried  in  the  direc¬ 
tion  parallel  to  the  plasma  flow.  Therefore,  for  simplicity,  we  omit 
effects  connected  with  electric  and  self-magnetic  fields.  Partially 
ionized  plasma  accelerates  in  the  axial  direction  due  to  the  pressure 
gradient  and  achieves  the  sound  speed  at  the  cavity  exit  plane.’ ^ 
Note  that  this  is  a  dominant  acceleration  mechanism  in  the  elec¬ 
trothermal  PPT-4  device,  whereas  traditional  PPT  has  predominant 
electromagnetic  acceleration  mechanism.  Teflon  evaporation  is  the 
origin  of  the  source  term  in  the  mass  conservation  equation.  In  this 
case,  the  governing  equations  in  dimensionless  form  are 


dv  P 

dz  «(1  —  v^) 


(1) 


dn  ndv  ^  p 

dz  vdz  V 


(2) 


where  p  =g/{NaCsma7trf^)  and  Cs  =  ([A:7e  -f* /: 7} ^  The  fol¬ 
lowing  boundary  conditions  are  used  for  Eqs.  (1)  and  (2):  n{z  = 
0)  =  1  and  du/dz  =0. 
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the  sheath.  Using  Maxwell’s  equation  and  assumption  3,  the  electric 
field  E{r)  changes  with  time  according  to  the  relation 

s^^  =  -[J.(r)  +  J,(r)]  (3) 

ot 

The  time  derivative  of  the  electric  field  is  a  function  of  the  radius  r. 
An  estimation  indicated^^  that  the  characteristic  time  for  the  electric 
field  to  reach  steady  state  decreases  with  radius  and  has  its  maximum 
value  at  the  MP  surface,  that  is,  where  r  =  Rp.  Thus,  Eq.  (3)  will  be 
solved  at  the  radius  Rp.  The  electron  and  ion  current  densities  are 
required  to  solve  Eq.  (3).  The  electron  current  density  is  calculated 
according  to  the  following  relation: 

4  =  J^exp(~eUp/kTe)  (4) 

The  ion  current  density  depends  on  the  ratio  of  Debye  length  to  the 
MP  radius  and,  thus,  will  be  different  for  small  and  large  MPs. 


Fig.  2  Plasma  density  and  velocity  spatial  distribution  in  the  cavity, 
13=0.6. 


The  plasma  density  and  velocity  distribution  are  shown  in  Fig.  2. 
Note  that  the  plasma  is  significantly  accelerated  toward  the  cavity 
exit.  Burton  et  al.  obtained  similar  flowfield  development  in  a 
liquid-injected  capillary  discharge^ ^  when  the  plasma  flow  ap¬ 
proaches  steady-state  conditions.  Taking  into  account  that  the 
plasma  velocity  should  be  sonic  at  the  exit  (z  =  1),  it  is  found  that 
p  —  0.6.  In  one  PPT-4  design  it  was  measured  that  the  ablation  rate 
is  about  30  /rg  per  lO-pts  pulse.^°  One  can  estimate  that  the  plasma 
density  near  the  central  electrode  should  be  m“^  using  this 

experimental  value  for  g.  Note,  however,  that  the  aforementioned 
estimated  number  density  of  the  plasma  in  the  cavity  represents  the 
low  limit.  This  is  because  it  was  assumed  that,  in  the  framework  of 
a  quasi-steady-state  model,  the  characteristic  time  is  equal  simply 
to  the  pulse  duration.  Under  the  assumption  that  the  characteristic 
discharge  time  is  about  the  acoustic  time  {L/C,,  which  is  typically 
a  few  microseconds)  the  plasma  density  near  the  central  electrode 
may  be  substantially  higher  (up  to  an  order  of  magnitude)  than  the 
preceding  estimation. 


B.  MP  Charging 

An  experimental  investigation  of  the  MP  size  distribution  showed 
a  variation  from  0.1  up  to  100  jLtm  (Ref.  5).  It  was  also  found  that 
some  MPs  have  a  spherical  shape,  whereas  generally  MPs  have  a 
rather  random  shape  distribution.^  To  make  it  possible  to  develop  a 
model  quantitatively,  we  will  assume  that  all  MPs  have  a  spherical 
shape  with  radii  in  the  experimentally  observed  range.  Under  the 
typical  conditions  of  PPT-4  operation,  the  electron  density  in  the 
cavity  during  the  pulse  is  about  lO^’-lO^"*  m“^,  which  corresponds  to 
a  debye  length  of  about  10“^- 10"’  m.  It  is  known  that  the  thickness 
of  the  sheath  around  a  MP  is  about  the  Debye  length.^* Therefore, 
one  can  see  that  there  are  basically  two  cases:  small  MPs  for  which 
the  Debye  length  is  less  than  or  about  the  MP  radius,  and  large  MPs 
for  which  the  Debye  length  is  much  less  than  the  MP  radius.  These 
two  cases  may  be  handled  using  different  approaches  for  modeling 
the  charging  process.  The  second  case  is  more  straightforward  and 
corresponds  to  that  of  a  plane  probe.  The  MP  charging  in  the  first 
case  was  developed  in  detail  in  Ref.  22. 

In  both  cases,  the  MP  charging  process  is  modeled  with  the  aid 
of  the  following  assumptions:  1)  The  plasma  consists  of  two  species 
of  charged  particles.  For  PPT-4  conditions,  thermal  equilibrium 
between  all  species  is  achieved  and,  therefore,  we  assume  that 
T  —  Te  =  Ti.  2)  The  plasma  Jet  flow  is  not  substantially  obstructed 
in  the  sheath  around  the  MP,  and,  thus,  spherical  symmetry  of  the 
plasma  density  relative  to  the  MP  is  assumed.  3)  Self-magnetic  field 
effects  are  neglected  because  PPT-4  has  predominantly  electrother¬ 
mal  acceleration  mechanism. 

The  kinetics  of  the  MP  charging  is  controlled  by  the  ion  and  elec¬ 
tron  fluxes  to  the  MP,  which  depend  on  the  potential  distribution  in 


].  Small  MPs 

This  case  corresponds  to  the  low  limit  of  the  measured  MP  size 
distribution  function.  In  the  case  of  a  spherical  sheath  around  an 
MP  and  L/,  >  the  ion  flux  may  be  calculated  using  the  orbital 
motion  limiF^  ^'*: 


Ji  =  Jiod-^eUp/kTi)  (5) 


This  expression  is  exact  for  an  infinite  Lp/Rp  ratio.  By  considering 
different  ion  trajectories  around  the  MP,  it  is  possible  to  calculate 
the  ion  flux  for  a  finite  Lp/Rp  ratio.  An  influence  of  this  effect  was 
examined  in  detail  in  Ref.  22. 

In  the  general  case,  the  capacitance  of  the  MP  in  the  plasma  is 
given  by  C  =  Qp/Up  =  47T RpSoG{Rp/Lp).  The  function  G(Rp/ 
L  p)  is  presented  in  Ref.  22,  and,  for  the  case  Rp/Lp  ^  1,  this  func¬ 
tion  is  G{Rp/Lp)^  1.8,  whereas  in  the  orbital  motion  limit  this 
function  is  equal  to  I. 

When  we  take  into  account  the  preceding  expression  for  MP 
capacitance  and  combine  Eqs.  (3-5),  the  kinetics  of  MP  charging 
may  be  described  by  the  following  dimensionless  equation  (in  the 
orbital  motion  limit): 


d0  Rp  1 


6 


dr  Lp 


+  (/>- 


(6) 


2.  Large  MPs 

In  this  case,  the  sheath  model  near  the  planar  electrode  can  be 
used.  Because  in  the  cavity  the  plasma  velocity  is  less  than  the 
sound  speed  (see  Sec.  II.A),  the  ion  flux  can  be  calculated  from  the 
Bohm  expression^^: 

4  =  OAeNC,  (7) 

The  capacitance  in  this  case  reads^^ 

C  =  47tRlso(\/R,  +  1/Ld)  ~  4nso{R'^/LD)  (8) 

When  we  take  into  account  the  expression  for  MP  capacitance,  the 
kinetics  of  MP  charging  may  be  described  by  the  following  dimen¬ 
sionless  equation: 


The  time  variation  of  the  dimensionless  MP  potential  is  shown  in 
Fig.  3  for  two  limiting  cases  of  Rp/Lp^^l  and  Rp/Lp  I .  Note 
that  the  charging  time  is  smaller,  whereas  the  steady-state  potential 
is  larger  in  the  case  Rp/Lp  1.  The  steady-state  potential  increases 
from  about  3.5  in  the  caseof  /?p/L£>  <5C  1  up  to  about  5  in  the  case  of 
Rp/Lp':^\.  All  possible  cases  realized  in  a  typical  PPT  are  placed 
between  the  limits  of  these  two  curves.  Thus,  in  the  PPT  plasma 
with  r  =  2  eV,  MP  has  a  negative  potential  of  about  -(7  -r  10)  V 
with  respect  to  the  plasma. 
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Fig.  3  Temporal  behavior  of  dimensionless  MP  potential  with  RplLo 
as  a  parameter. 


carrying  plasma.  To  predict  the  MP  velocity,  we  simply  integrate  the 
equation  of  motion  from  some  starting  point.  We  assumed  that  MPs 
have  zero  initial  velocity.  Given  the  solution  for  plasma  density  and 
velocity  distribution  in  the  cavity,  we  evaluate  the  forces  that  act  on 
an  individual  MP.  The  equation  of  motion  of  an  individual  MP  may 
be  written  as 


mpVp) 

dt 


F’d  +  QpE 


(11) 


The  mass  Mp  of  the  MP  traveling  in  the  plasma  may  be  changed  as 
a  result  of  MP  ablation  (see  next  section). 

The  first  term  is  the  total  drag  force  and  the  second  term  is  the 
electric  force,  which  depends  on  the  electric  field  in  the  plasma  and 
the  MP  charge.  An  average  electric  field  E  in  the  plasma  may  be  cal¬ 
culated  from  Ohm’s  law  for  known  current  density  7,  namely,  E  = 
j/a.  We  will  consider  a  diffuse  discharge  on  the  anode  and,  there¬ 
fore,  j  =  I / TV The  current  density  was  estimated  at  the  current 
peak  of  about  8  kA  and  plasma  conductivity  was  calculated  for 
electron  temperature  of  about  2  eV  (see  Ref.  20). 

In  partially  ionized  plasmas,  there  are  coulomb  collisions  of  the 
MP  with  ions  because  the  MP  is  charged.  These  collisions  con¬ 
tribute  to  the  drag  and  may  result  in  a  force  called  ion  drag.^^  In 
low-density  plasmas  with  large  ratio  of  Debye  length  to  MP  radius 
(small  MPs,  see  preceding  section),  it  was  found  that  the  momen¬ 
tum  transfer  due  to  ions  that  are  collected  by  a  MP  cross  section 
is  much  less  than  that  due  to  ions  scattered  but  not  collected.^* 
Note  that  if  the  Debye  length  is  much  less  than  or  about  the  MP 
size  (large  MPs),  the  momentum  transfer  is  determined  by  the  MP 
cross  section.^^  In  general,  the  neutral  drag  force  is  also  determined 
by  the  MP  cross  section.  Analyses  of  the  MP  size  distribution  and 
range  of  possible  plasma  densities  (10^^ -lO^"*  m“^)  shows  that  both 
free  molecular  and  continuum  regimes  for  the  drag  force  can  be 
realized. 


7.  Small  MP  ( Free  Molecular  Regime } 

In  the  free  molecular  regime,  that  is,  7?^  <$C  the  drag  force  can 
be  written  as^^ 

Fj  =  {^pRl /2^)[[s  +  (l/2i)  expC-i^)] 


Fig.  4  MP  char;^  dependence  on  the  radius  with  electron  density  as  a 
parameter,  =  1.5  eV. 

5.  MP  Charge 

To  calculate  possible  values  of  the  charge  accumulated  on  the  MP, 
we  will  convert  from  the  dimensionless  parameters. The  MP  charge 
(in  the  large  MP  limit)  may  be  connected  with  the  dimensionless 
potential  as  follows: 


-i-(5^  +  1  -  l/4.y^)x/7rerf(5)}  (12) 

where  p  =ml2kT  and  s  —  2Vpl[-yJ{Tt){V  —  V^)]. 

2.  Large  MP  (Continuum  Regime) 

In  this  case  the  drag  force  that  acts  on  the  isolated  MP  placed  into  a 
plasma  moving  with  velocity  V  is  determined  by  the  expression^^’^* 


Q,  =  4nRlLDN4>  (10a) 

In  the  small  MP  limit,  the  corresponding  expression  for  the  MP 
charge  reads 

Qp  =  47r  SoRp(f>kTJe  (10b) 

The  calculated  steady-state  MP  charge  as  a  function  of  MP  radius 
is  shown  in  Fig.  4  with  electron  density  as  a  parameter  for  the  large 
MP  limit  (7?^  >  1  /um).  Note  that  the  Ujim  MP  has  a  charge  of  about 
10^  electrons  in  the  case  of  N  =  10^^  m"^.  One  can  see  that,  in  this 
case,  MP  charge  increases  with  electron  density  unlike  the  other 
case  (small  MP  with  Rp  <  1  fim)  when  Qp  is  independent  of  the 
plasma  density.  In  the  latter  case,  the  MP  charge  is  a  linear  function 
of  Rp  and  can  be  estimated  as  Qp  =  2.4  •  10^  x  7?;,  (micrometers) 
where  Qp  is  in  electrons. 

C.  Momentum  Transfer 

There  are  two  primary  forces  acting  on  the  MP  in  the  PPT  dis¬ 
charge:  1)  a  drag  force  due  to  collision  with  neutral  atoms  and  ions 
and  2)  a  force  due  to  the  presence  of  the  electric  field  in  a  current 


F,  =  CD7rRlpiV-V,ff2  (13) 

where  C/)  is  the  drag  coefficient  that  depends  on  the  Reynolds  num¬ 
ber.  For  instance,  in  the  case  of  Re  <  1,  coefficient  Co  =24.8/7?e 
(Refs.  30  and  31).  Estimation  shows  that  in  the  PPT  plasma,  the 
Reynolds  number  is  about  10“^  V,  where  V  is  the  plasma  velocity 
in  meters  per  second. 

The  calculation  of  the  MP  velocity  in  the  continuum  regime  as  a 
function  of  time  with  normalized  plasma  velocity  as  a  parameter  is 
plotted  in  Fig.  5a  for  the  case  Rp  =  1  ^m.  Note  that  the  maximum 
possible  MP  velocity  is  about  230  m/s.  In  experiments  it  was  found 
that  some  MPs  have  a  velocity  of  about  200  m/s  (Ref.  5).  Smaller 
MPs  may  have  a  larger  velocity  as  shown  in  Fig.  5b  for  the  case  of 
Rp—O.l  /xm,  calculated  in  the  free  molecular  regime.  One  can  see 
that  the  MP  velocity  depends  on  the  residual  time  in  the  discharge. 
Thus,  those  MPs  generated  toward  the  end  of  the  pulse  are  expected 
to  have  a  smaller  velocity. 

The  ratio  of  MP  velocity  to  the  plasma  velocity  dependence  on 
the  MP  radius  is  shown  in  Fig.  6  with  plasma  density  as  a  parameter. 
These  calculationsare  performed  for  the  case  of  MP  residual  time 
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D.  Heat  Transfer 

As  the  particle  moves  through  the  plasma,  it  is  heated  by  neutral, 
electron,  and  ion  fluxes  and  is  cooled  by  radiation.  An  additional 
cooling  mechanism  is  due  to  decomposition  of  the  dielectric  material 
under  high  temperature.  We  will  consider  the  situation  where  the 
thermal  time  constant  is  much  shorter  than  the  MP  residual  time.  In 
this  case,  the  MP  has  a  uniform  temperature  throughout  its  volume 
and  inward  heat  conduction  from  the  surface  may  be  neglected.  The 
resulting  energy  balance  reads 

4  dT 

-ttRIpC  +  Q„-Q,-  Q,)  (14) 

According  to  the  charging  model  (see  Sec.  II .B)  the  ion  and  electron 
fluxes  are  equal  afterthe  MP  reaches  steady-state  charge.  Therefore, 
the  gross  energy  input  by  ions  and  electrons  to  the  MP  is  given  by  the 
current  density  and  the  sum  of  the  energy  carried  by  each  species: 


h)Rp=0Afim  time,  jus 

Fig.  5  Temporal  behavior  of  the  MP  velocity  with  dimensionless 
plasma  velocity  as  a  parameter,  Tg  =  1.5  eV  and  N  =  10^^  m”^. 


Particle  size,  Rp  pm 


Fig.  6  Normalized  MP  velocity  dependence  on  the  MP  radius  with 
plasma  density  as  a  parameter,  Tg  =  1.5  eV  and  At  =  10  fis. 


Qi,e  =  j,{Rp)[2kTe+kTg-^eUp{t)]  (15) 

where  Up(t)  is  the  time-dependent  MP  potential  (see  Sec.  II. A). 

The  neutral  atoms’  heat  transfer  may  be  calculated  using  the  New¬ 
tonian  model 

Qn^^cxAT  (16) 


In  the  general  case,  the  coefficient  of  heat  transfer  is  a  complex 
function  of  MP  size,  gas  flow  temperature  and  velocity,  heat  conduc¬ 
tivity,  specific  heat,  and  density.  The  relation  between  the  coefficient 
of  heat  transfer  and  the  preceding  parameters  was  determined  by  a 
similarity  law.^°  For  the  convective  heat  transfer  between  a  body 
and  gas  flow,  the  following  similarity  can  be  used: 

Nu  =  aRp/k  (17) 

Considering  plasma- MP  interaction  in  the  PPT-4  cavity,  one  can 
estimate  Nusselt  number  in  the  case  when  the  particle  is  not  moving 
relative  to  the  plasma.  In  this  case  the  Nusselt  number  Nu  equals  2 
(Ref.  30).  The  dependence  of  atomic  thermal  conductivity  on  tem¬ 
perature  is  given  by:  A  =  2.4  x  W/mK  (Ref.  30). 

The  radiative  flux  is  given  by  the  Stefan- Boltzmann  law: 

Qr=<^r^T^  (18) 

The  heat  flux  associated  with  material  decomposition  can  be  cal¬ 
culated  as 

Q,  =  V,AH  (19) 

The  ablation  rate  can  be  estimated  at  equilibrium  using  Knudsen’s 
law: 


=  Po/P 


/  niaf  (IJxkT) 


(20) 


where  Pq  is  the  equilibrium  pressure  of  Teflon  (see  Refs.  1  and  4). 

The  rate  of  the  MP  radius  change  is  determined  by  the  ablation 
rate  : 


(21) 


equal  to  the  discharge  pulse,  that  is,  for  /  =  10  fis.  One  can  see  that 
small  MPs  in  the  dense  plasma  can  approach  the  plasma  velocity 
and  can  be  entrained  by  the  plasma  jet.  However,  the  large  MPs  are 
generally  slower  and  flow  substantially  behind  the  plasma  jet. 

Note  that,  due  to  ablation,  the  MP  radius  decreases  in  the  course 
of  interaction  with  the  plasma,  and,  thus,  MP-plasma  interaction 
may  change  from  initially  continuum  to  free  molecular.  An  effect 
of  MP  ablation  is  considered  in  the  next  section. 


The  initial  condition  is  Rp{t  =  0)  =  Rp^. 

The  calculated  temporal  variation  of  the  MP  temperature  is  shown 
in  Fig.  7  with  MP  radius  as  a  parameter.  It  can  be  seen  that  small  MPs 
are  heated  substantially  up  to  1000  K  during  a  short  time  period, 
whereas  large  MPs  are  only  slightly  affected  by  the  plasma.  Heating 
of  the  MP  leads  to  significant  ablation  as  plotted  in  Fig.  8.  Note  that 
small  MPs  (~  1  fjim)  completely  ablate  if  flieir  residual  time  is  larger 
than  0.5  /as. 
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Fig.  7  Temporal  dependence  of  the  MP  temperature  with  the  MP 
radius  as  a  parameter,  7,  =  1.5  eV  and  N  =  10^^  m“^. 


Fig.  8  MP  radius  vs  time  with  the  initial  MP  radius  as  a  parameter, 
7,  =1.5  eV  and  A^  =  10^3 


III.  Discussion 

One  of  the  important  possible  applications  of  the  PPT  is  for  at¬ 
titude  control  of  small  satellites  that  are  part  of  a  group  of  close¬ 
flying  satellites.  In  this  case,  interaction  of  the  integrated  plasma 
plume  from  the  PPT  with  another  spacecraft  becomes  an  impor¬ 
tant  issue.  An  experimental  investigation  has  shown  that  the  plasma 
plume  contains  a  substantial  fraction  of  MPs  that  originated  from 
the  dielectric  as  well  as  from  the  electrodes.  The  presence  of  such 
particles  in  the  interacting  plume  may  be  a  source  of  local  damage 
to  the  spacecraft.  The  potential  for  MP  damage  on  their  collision 
with  spacecraft  surfaces  depends  on  the  MP  directed  velocity.  The 
present  calculation  shows  that  the  MP  velocity  depends  mainly  on 
MP  size  and  residual  time  in  the  discharge.  It  was  assumed  that 
MP- plasma  interaction  occurs  mainly  during  the  discharge  pulse. 
This  is  because  the  plasma  density  and  temperature  that  control  the 
interaction  drop  substantially  after  the  pulse.  Therefore,  the  maxi¬ 
mum  MP  residual  time  considered  in  this  model  was  equal  to  the 
pulse  duration.  Those  small  MPs  emitted  at  the  beginning  of  the 
discharge  pulse  may  have  a  velocity  of  200  m/s.  Thus,  one  can  ex¬ 
pect  a  relatively  broad  range  of  MP  velocity  from  zero  up  to  a  few 
hundred  meters  per  second  depending  on  the  point  of  MP  forma¬ 
tion  and  the  MP  size.  Another  useful  outcome  from  the  MP  velocity 
calculation  is  the  possibility  to  estimate  some  characteristic  time  of 
MP  flux  ejection.  Because  MP  velocity  is  generally  smaller  than 
the  plasma  velocity,  and  due  to  the  broad  range  of  MP  velocities, 
the  effective  time  of  MP  ejection  from  the  thruster  is  much  larger 
than  the  pulse  duration.  For  instance,  our  estimation  shows  that  the 
large  MPs  (~100  fxm)  exhaust  from  the  thruster  only  after  several 
seconds. 


Because  of  the  possible  electrostatic  nature  of  MP  interaction 
with  a  spacecraft,  the  MP  charge  may  play  an  important  role.  A 
spacecraft  interacts  with  a  plasma  jet  so  as  to  become  electrically 
charged  to  ensure  zero  net  current,  just  like  the  MP.  Thus,  because  the 
spacecraft  has  some  negative  potential  with  respect  to  the  plasma, 
the  MP  may  interact  with  the  electric  field  in  the  thin  sheath  around 
the  spacecraft.  In  some  cases,  the  MP  may  be  reflected  in  the  sheath 
without  mechanical  collision  with  the  surface.  This  effect  has  a 
positive  role  because  the  electrostatic  nature  of  MP  interaction  with 
surfaces  in  space  is  not  then  connected  with  any  mechanical  damage. 
Thus,  it  is  important  to  estimate  the  possibility  of  such  MP- surface 
interactions.  Previously,  this  effect  was  studied  for  the  case  of  vac¬ 
uum  arc-generated  MPs.  It  was  obtained  both  theoretically  and 
experimentally  for  typical  conditions  that  the  electrostatic  nature  of 
the  MP  interaction  with  surfaces  (without  mechanical  touching)  is 
possible.  It  occurs  in  those  cases  when  an  MP  approaches  the  surface 
with  a  grazing  angle  such  that  the  MP  kinetic  energy  will  be  smaller 
than  the  potential  electrostatic  energy.  The  possibility  of  this  effect 
may  be  estimated  using  the  following  integral  expression: 

QpU.  >  M,V^j2 

For  instance,  if  the  spacecraft  potential  Us  is  about  10  V,  a  l-jxm  MP 
may  be  electrostatically  reflected  if  its  normal  velocity  component 
is  less  than  4  m/s. 

Another  important  outcome  from  the  present  study  is  prediction 
of  the  MP  complete  or  partial  decomposition  due  to  interaction  with 
the  plasma  during  the  discharge  pulse.  It  was  found  that  MPs  with 
size  less  than  10  /xm  may  totally  decompose  during  10  jjls.  However, 
large  MPs  may  be  still  be  present  in  the  plasma  plume.  This  means 
that  the  size  distribution  of  MPs  exhausted  from  the  PPT  is  different 
from  the  original  MP  distribution.  According  to  our  calculation,  the 
difference  is  more  significant  in  the  range  of  small  MPs.  Thus,  any 
measurements  of  the  MP  flux  in  the  plasma  plume  involve  some 
combination  of  the  original  MP  distribution  and  that  from  MP- 
plasma  interaction  rather  than  the  size  distribution  of  MPs  emitted 
from  the  Teflon.^^  Note  also  that  the  MP  decomposition  may  be  a 
substantial  volume  source  of  the  neutral  component  of  the  discharge 
plasma.  However,  because  there  is  an  uncertainty  in  the  original 
MP  size  distribution,  the  significance  of  the  last  effect  cannot  be 
estimated  accurately.  One  can  also  conclude  that  the  increase  in  MP 
residence  time  leads  to  MP  decomposition  and,  therefore,  may  have 
some  effect  on  increasing  propellant  efficiency. 

Note  that  in  this  paper  the  MP-plasma  interaction  was  considered 
for  the  conditions  typical  for  the  electrothermal  PPT-4  device.  This 
device  is  fundamentally  different  from  the  usual  LES  8/9  class  PPT 
with  dominant  electromagnetic  acceleration  mechanism.  It  appears 
that  the  basic  features  of  the  present  model  should  be  same  for 
electromagnetic  PPT.  The  main  difference  is  a  much  higher  plasma 
velocity  realized  in  this  PPT  due  to  J  x  B  acceleration.  It  implies 
that  the  effective  MP  residual  time  may  decreases,  and,  therefore, 
the  MP  charging,  acceleration,  and  heating  may  be  affected.  Higher 
plasma  velocity  will  increase  the  ion  flux  to  the  MP  that  may  lead 
to  increase  MP  charge  in  the  steady  state.  However,  because  the 
MP  steady-state  charge  depends  logarithmically  on  the  ion  flux, 
one  would  expect  only  small  changes.  Small  (submicrometer)  MPs 
have  the  potential  for  acceleration  when  they  interact  with  a  high- 
velocity  plasma  jet. 

IV.  Conclusion 

We  have  shown  that  the  plasma  may  affect  particulate  contami¬ 
nation  from  the  plume  of  a  PPT.  During  flight,  a  MP  may  become 
charged,  heated,  and  accelerated. It  was  found  that  the  charging  time 
depends  on  MP  size  and  is  generally  smaller  in  the  case  of  small 
MPs.  We  have  found  that  a  l-/xm  MP  has  a  charge  of  about  10^ 
electrons  in  the  case  of  an  electron  density  of  about  10^^  m“^.  There 
are  two  primary  forces  acting  on  the  MP  in  the  PPT  discharge:  a  drag 
force  due  to  collision  with  neutral  atoms  and  ions  and  an  electric 
force  due  to  the  presence  of  the  electric  field  in  a  current  carry¬ 
ing  plasma.  The  maximum  possible  velocity  of  a  l-;xm  MP  is  about 
230  m/s,  which  is  close  to  that  estimated  experimentally.  Only  small 
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MPs  (~0.1  fjim)  can  be  entrained  by  the  plasma  jet,  whereas  large 
MPs  are  generally  slower  and  flow  substantially  behind  the  plasma 
jet.  Small  MPs  (<  1  fxm)  may  completely  decompose  during  a  l-jis 
pulse. 
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Abstract 

In  this  work  we  present  a  model  of  the  near  field  plasma  plume  of  a  Pulsed  Plasma  Thruster  (PPT).  As  a 
working  example  we  consider  a  micro-PPT  developed  at  the  Air  Force  Research  Laboratory.  This  is  a 
miniaturized  design  of  the  axisymmetric  PPT  with  a  thrust  in  the  10  pN  range  that  utilizes  Teflon™  as  a 
propellant.  The  plasma  plume  is  simulated  using  hybrid  fluid-PIC-DSMC  approach.  The  plasma  plume 
model  is  combined  with  Teflon  ablation  and  plasma  generation  models  that  provide  boundary  conditions 
for  the  plume.  This  approach  provides  a  consistent  description  of  the  plasma  flow  from  the  surface  into  the 
near  plume.  The  magnetic  field  diffusion  into  the  plume  region  is  also  considered  and  plasma  acceleration 
by  the  electromagnetic  mechanism  is  studied.  Teflon  ablation  and  plasma  generation  analyses  show  that 
the  Teflon  surface  temperature  and  plasma  parameters  are  strongly  non-uniform  in  the  radial  direction.  The 
plasma  density  near  the  propellant  surface  peaks  at  about  10^"^  m"^  in  the  middle  of  the  propellant  face  while 
the  electron  temperature  peaks  at  about  4.5  eV  near  the  electrodes.  The  model  predicts  ablated  mass  per 
pulse  of  about  1  \x,g  that  is  close  to  that  measured  in  experiment.  The  plume  simulation  shows  that  a  dense 
plasma  focus  is  developed  at  a  few  mm  from  the  thruster  exit  plane  at  the  axis.  This  plasma  focus  exists 
during  the  entire  pulse,  but  the  plasma  density  in  the  focus  decreases  from  about  2x10^^  m"^  at  the 
beginning  of  the  pulse  down  to  0.3x10^^  m"^  at  5  ps.  The  velocity  phase  is  centered  at  about  30  km/s  in  the 
axial  direction.  At  later  stages  of  the  pulse  there  are  two  ion  populations  with  positive  and  negative  radial 
velocity.  An  ion  population  having  negative  axial  velocity  up  to  about  10  km/s  is  predicted.  This  is  a 
significant  finding  in  terms  of  backflow  contamination  onto  a  spacecraft. 
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1.  Introduction 


The  pulsed  plasma  thruster  (PPT)  was  among  the  first  of  various  electrical  propulsion  concepts 
accepted  for  space  flight  mainly  due  to  its  simplicity  and  hence  high  reliability'.  However,  the 
PPT  has  an  efficiency  that  is  generally  low^  at  about  10%  leaving  open  the  opportunity  for 
considerable  improvement^.  Currently,  PPT’s  are  considered  as  an  attractive  propulsion  option 
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for  stationkeeping  and  drag  makeup  purposes  of  mass  and  power  limited  satellites  ’ . 
Guaranteeing  successful  operation  of  spacecraft  using  a  PPT  requires  a  complete  assessment  of 
the  spacecraft  integration  effects.  The  PPT  plume  contains  various  ion  and  neutral  species  due  to 
propellant  decomposition  and  possible  electrode  erosion.  Some  attempts  of  PPT  plume  modeling 
using  particle  simulations  were  performed  recently*’’’*.  In  Refs.  7,8  we  have  considered  the  plume 
flowfield  exhaust  from  an  electrothermal  PPT  and  therefore  electromagnetic  effects  in  the  plume 
were  neglected.  Different  variations  of  electromagnetic  PPTs  are  also  candidates  for  various 
missions’.  Recently,  a  micro-PPT  has  being  designed  at  AFRL  for  delivery  of  very  small  impulse 
bit'°.  This  is  a  simplified  miniaturized  version  of  a  conventional  PPT  designed  to  provide  attitude 
control  and  stationkeeping  for  microsatellites.  We  will  use  the  AFRL  micro-PPT  as  a  working 
example  for  several  reasons.  Firstly,  electromagnetic  OxB)  acceleration  is  the  primary 
mechanism  in  this  thruster;  and  secondly,  there  is  no  internal  flow  in  this  device  and  therefore  the 
near-field  plasma  plume  is  an  essential  part  of  the  thrust  generation  process.  Therefore  careful 
modeling  of  the  acceleration  is  needed  to  understand  the  characteristics  of  the  device  as  a  whole 
in  addition  to  being  a  pre-cursor  to  accurate  estimation  of  contamination  issues.  Since  in  this 
device  there  is  no  separation  between  the  main  plasma  acceleration  region  and  the  plume 
expansion,  both  regions  must  be  simulated  in  one  model.  Because  the  plasma  acceleration  is 
external,  the  plasma  is  sufficiently  rarefied  so  that  an  MHD  approach  such  as  MACH2  (Ref.  1 1) 
cannot  be  used. 
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An  accurate  model  of  the  PPT  plume  relies  on  the  boundary  and  initial  conditions.  These 
conditions  can  be  formulated  by  consideration  of  the  Teflon  ablation  process.  The  Teflon  ablation 
computation  is  based  on  a  recently  developed  kinetic  ablation  model'^’'^.  In  this  model  the 
detailed  physics  of  the  Teflon  evaporation  is  studied  by  consideration  of  the  distribution  function 
of  the  particles  in  the  kinetic  layer  adjacent  to  the  surface. 

Another  important  effect  related  to  the  plasma  plume  exhaust  from  an  electromagnetic  PPT  is  the 
magnetic  field  diffusion  into  the  near  plume.  Previously,  we  have  modeled  the  effect  of  the 
magnetic  field  on  the  near-field  plume  for  Hall  thrusters'"*  under  steady  state  conditions.  It  was 
found  that  the  magnitude  of  the  magnetic  field  at  the  thruster  exit  has  an  important  effect  on  the 
plasma  potential  distribution  in  the  plume.  In  the  present  research,  it  is  proposed  to  include  the 
electromagnetic  effects  on  the  near  field  plume  of  unsteady  plasma  flow.  The  computational 
domain  is  shown  in  Fig.  1.  The  model  is  based  on  a  hybrid  approach  involving  a  DSMC 
description  of  neutrals,  a  PIC  model  for  ions  and  a  fluid  description  of  the  electrons.  In  these 
methods,  the  potential  distribution  is  usually  calculated  by  reducing  the  electron  momentum 
equation  to  the  Boltzmann  relation  in  the  absence  of  a  magnetic  field.  In  the  plasma  plume 
domain  where  the  magnetic  field  exists,  i.e.  the  near  field,  it  is  necessary  to  include  the  magnetic 
field  effects  in  the  electron  momentum  equation. 

2.  The  model  of  the  plasma  layer 

The  model  presented  here  describes  the  plasma  layer  near  the  Teflon  surface  as  shown  in  Fig.  2. 
The  model  of  the  plasma  layer  includes  Joule  heating  of  the  plasma,  heat  transfer  to  the  Teflon, 
and  Teflon  ablation.  Mechanisms  of  energy  transfer  from  the  plasma  column  to  the  wall  of  the 
Teflon  include  heat  transfer  by  particle  convection  and  by  radiation.  The  Teflon  ablation 
computation  is  based  on  a  recently  developed  kinetic  ablation  model'^.  It  is  assumed  that  within 
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the  plasma  layer  all  parameters  vary  in  the  radial  direction  r  (see  Fig.  2).  The  energy  balance 
equation  can  be  written  in  the  form: 

|nedTe/dt  =  Qj-Qr-Qp . (1) 

where  Qj  is  the  Joule  heat,  Qr  is  the  radiation  heat  and  Qp  is  the  heat  associated  with  particle 
fluxes.  This  equation  depends  on  the  coordinate  along  the  propellant  face.  For  known  plasma 
density  and  temperature  the  heat  flux  to  the  surface  is  calculated.  The  Teflon  surface  temperature 
is  calculated  from  the  heat  transfer  equation  with  boundary  conditions  that  take  into  account 
vaporization  heat  and  conductivity.  The  solution  of  this  equation  is  considered  for  two  limiting 
cases  of  substantial  and  small  ablation  rate  very  similar  to  that  described  in  Ref.  8.  The  density  at 
the  Teflon  surface  is  calculated  using  the  equilibrium  pressiue  for  Teflon.  The  plasma  density  in 
the  layer  is  determined  in  the  framework  of  the  kinetic  ablation  model  (see  next  section).  For 
known  pressure  and  electron  temperature  one  can  calculate  the  chemical  plasma  composition 
assuming  LTE*’'^’’*.  The  Saha  equations  are  supplemented  by  the  conservation  of  nuclei  and 
quasi-neutrality. 

3.  Ablation  model 

The  Teflon  ablation  is  modeled  in  the  framework  of  the  approximation'^  based  on  a  kinetic  model 
of  the  material  evaporation  into  discharge  plasmas'^.  The  model  couples  two  different  layers 
between  the  surface  and  the  plasma  bulk  as  shown  in  Fig.  2b:  (1)  a  kinetic  non-equilibrium  layer 
adjusted  to  the  surface  with  a  thickness  of  about  one  mean  free  path;  and  (2)  a  collision- 
dominated  layer  with  thermal  and  ionization  non-equilibrium.  The  velocity  at  the  edge  of  the 
kinetic  layer  U|  can  be  determined  from  the  coupling  solution  of  the  hydrodynamic  layer  and  the 
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quasi-neutral  plasma.  For  known  velocity  and  density  at  this  interface,  it  is  possible  to  calculate 
the  ablation  rate.  In  the  hydrodynamic  layer  the  relation  between  the  velocities,  temperatures  and 
densities  at  the  boundaries  1  and  2  as  well  as  the  ablation  rate  are  formulated  according  to  Ref.  13 
in  the  form: 


r  =  mU|N,  =  N,[(2kT,/m)-(T2N2/2T,  -N,/2)/(N,-N,"/N2)f  ^ . (2) 

The  system  of  equations  is  closed  if  the  equilibrium  vapor  pressure  can  be  specified  that 
determines  parameters  (No  and  To)  at  the  Teflon  surface.  The  full  self-consistent  solution  of  this 
problem  can  be  obtained  when  the  ablation  is  coupled  with  the  plasma  plume  expansion.  In  the 
present  work  in  order  to  simplify  the  problem,  we  will  assume  that  the  plasma  accelerates  up  to 
the  sound  speed  near  the  boundary  2.  This  assumption  can  be  justified  by  the  fact  that  due  to 
significant  electrodynamic  acceleration  in  this  type  of  PPT,  the  plasma  density  will  quickly 
decrease  therefore  providing  solution  of  the  ablation  problem  close  to  that  ablation  into  the 
vacuum.  In  this  case  the  plasma  density  at  the  edge  of  the  kinetic  layer  will  be  equal  to  0.34  No 
and  the  temperature  is  0.7  To.  The  flux  returned  to  the  surface  is  equal  to  16%  of  the  ablated  flux 
(Ref  12). 

4.  Plasma  plume  electrodynamics 

The  general  approach  for  the  plume  model  is  based  on  a  hybrid  fluid-particle  approach  that  was 
used  previously  (Refs.  7).  In  this  model,  the  neutrals  and  ions  are  modeled  as  particles  while 
electrons  are  treated  as  a  fluid.  Elastic  (momentum  transfer)  and  non-elastic  (charge  exchange) 
collisions  are  included  in  the  model.  The  grids  employed  in  this  computation  are  also  similar  to 
those  used  previously  (Ref. 7).  The  particle  collisions  are  calculated  using  the  direct  simulation 
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Monte  Carlo  (DSMC)  method’’.  Momentum  exchange  cross  sections  use  the  model  of  Dalgamo 
et  al.’*,  while  charge  exchange  processes  use  the  cross  sections  proposed  by  Sakabe  and  Izawa'®. 
Acceleration  of  the  charged  particles  is  computed  using  the  Particle-In-Cell  method  (PIC)’“.  The 
plasma  velocity  distribution  depends  upon  the  magnetic  field  distribution  and  ion  dynamics  is 
calculated  as  follows: 

dV/dr  =  -Cs’Vln(n)  +  jxB/mn . (3) 

where  Cs  is  the  soimd  speed,  n  is  the  plasma  density,  j  is  the  current  density  and  B  is  the  magnetic 
field. 

The  electron  dynamics  is  very  important  in  the  plasma  plume.  Previously  our  model  was  based  on 
the  assumption  that  electrons  rapidly  reach  the  equilibrium  distribution  and  in  the  absence  of  the 
magnetic  field  can  be  described  according  to  the  Boltzmann  distribution.  While  this  was  a 
satisfactory  assumption  in  the  case  of  an  electrothermal  thruster  plume  this  is  not  suitable  for  the 
near  field  of  an  electromagnetic  thruster.  In  the  presence  of  a  strong  magnetic  field,  the  electron 
density  distribution  deviates  from  that  according  to  Boltzmann”.  In  the  case  of  a  magnetic  field 
the  electron  momentum  equation  reads  (neglecting  electron  inertia): 

0  =  -e’ne(E+VeXB)  -  eVPg  -  Veimj . (4) 

We  have  assumed  quasi-neutrality  therefore  ne  =  nj  =  n.  The  electric  and  magnetic  field 
distributions  in  the  plume  are  calculated  from  the  set  of  Maxwell  equations.  We  further  assume 
that  the  magnetic  field  has  only  an  azimuthal  component  and  also  neglect  the  displacement 
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current.  The  combination  of  the  Maxwell  equations  and  electron  momentum  conservation  gives 
the  following  equation  for  the  magnetic  field: 


dB/at  =  1/(cth)V^B  -  Vx(jxB/(en))  +  Vx(VxB) 


(5) 


where  a  is  the  plasma  conductivity,  |X  is  the  permittivity,  n  is  the  plasma  density,  j  is  the  current 
density  and  V  is  the  plasma  velocity.  A  scaling  analysis  shows  that  the  various  terms  on  the  right 
hand  side  of  Eq.  5  may  have  importance  in  different  regions  of  the  plasma  plume  and  therefore 
general  end-to-end  plasma  plume  analysis  requires  keeping  all  terms  in  the  equation.  In  the  case 
of  the  near  plume  of  the  micro-PPT  with  a  characteristic  scale  length  of  about  1  cm  the  magnetic 
Reynolds  number  Rem«l  and  therefore  the  last  term  can  be  neglected.  Taking  this  into  account 
in  the  dimensionless  form,  Eq.  5  can  be  written  as: 


at 


Re„,  — =V"B 


•(0JT>{ 


a  Br  a(Br) 
arV  n^  dz 


a  Br^ 

azT-m  dr  ^ 


(6) 


where  (cox)  is  the  Hall  parameter  that  measures  of  the  Hall  effect.  Therefore,  depending  on  the 
plasma  density,  the  Hall  effect  may  be  important  for  the  magnetic  field  evolution.  One  of  the  first 
calculations  of  the  plasma  flow  with  Hall  effect  were  performed  by  Brushlinski  and  Morozov  (see 
Ref.  22  and  references  therein).  They  considered  isothermal  flow.  The  plasma  density  becomes 
high  at  the  cathode  and  lower  at  the  anode.  The  Hall  effect  has  a  particularly  noticeable  influence 
on  the  magnetic  field  distribution.  The  field  near  the  anode  increases  and  near  the  cathode 
decreases.  As  a  result  the  current  is  deflected  to  the  side  and  grazes  the  anode. 

Our  estimations  show  that  the  Hall  parameter  ci>T«l  if  the  plasma  density  near  the  Teflon  surface 
N>10^^  m'^.  This  case  is  realized  in  the  micro-PPT  (see  the  next  section)  so  the  Hall  effect  is 
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expected  to  be  small  for  this  particular  case.  Having  the  magnetic  field  distribution  one  can 
calculate  the  current  density  distribution  from  Ampere’s  law: 

lij  =  VxB . (7) 

The  magnetic  field  and  current  distributions  calculated  from  this  model  are  used  in  PIC  to 
evaluate  the  ion  dynamics  according  to  Eq.  3. 

5.  Boundary  conditions 

The  boundary  conditions  for  the  magnetic  field  calculations  are  shown  in  Fig.  1.  We  have 
assumed  that  the  current  is  uniform  on  both  electrodes  that  allows  us  to  estimate  the  current 
density  on  the  cathode  jc  and  on  the  anode  ja.  The  magnetic  field  is  assumed  to  vary  as  1/r  on  the 
upstream  boundary.  At  the  lateral  boundary  we  have  assumed  that  the  normal  current  jn=0.  The 
downstream  boundary  is  considered  to  be  far  enough  away  that  B=0  can  be  assumed.  Along  the 
centerline  the  magnetic  field  B=0. 

The  boundary  conditions  for  the  plume  are  generated  through  solution  of  the  Teflon  ablation 
problem  as  will  be  presented  in  the  Results  section.  These  are  time  and  radial  dependent 
variations  of  the  plasma  (including  Carbon  and  Flourine  ions  and  neutrals)  density  and  electron 
temperature. 

The  results  are  presented  for  a  3.6  mm  diameter  micro-PPT  which  has  a  0.9  mm  diameter  central 
electrode,  3.1  mm  propellant  diameter  and  0.24  mm  anode  wall  (Ref  10).  In  these  simulations. 
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the  experimental  current  waveform  was  used,  that  is  described  in  a  first  approximation  as  an 
underdamped  LRC  circuit  current: 


I(t)  =Ip-sin(at)exp(-Pt) 


where  Ip  = 


P=  L  is  the  effective  inductance  in  the  circuit,  C  is 


the  capacitance,  R  is  the  total  circuit  resistance,  and  E  is  the  pulse  energy.  The  best  fit  with  the 
experimental  waveform  (frequency)  corresponds  to  a=3-10*  s  '.  For  C=0.3  pF  we  can  estimate 
that  L  in  the  circuit  is  about  3.6- lO  ’  H.  Results  presented  below  correspond  to  the  15.2  J  (Ref  10). 


6.  Results 

The  spatial  and  temporal  variation  of  the  Teflon  surface  temperature  is  shown  in  Fig.  3a.  The 
Teflon  temperature  sharply  increases  during  the  first  1-2  ps  of  the  pulse  and  peaks  at  about  635 
K.  One  can  see  that  the  temperature  is  generally  non-uniform  in  the  radial  direction  and  has  a 
minimum  at  radial  distances  of  1.1 -1.3  mm.  Since  the  Teflon  ablation  is  approximately 
exponentially  proportional  to  the  surface  temperature,  the  model  predicts  a  lower  rate  of  ablation 
in  the  areas  where  the  surface  temperature  has  a  minimum.  Taking  this  into  account,  the  effect  of 
the  temperature  distribution  may  be  related  to  the  preferential  charring  of  the  Teflon  siuface 
observed  experimentally  [Ref.  10].  A  detailed  study  of  the  Teflon  surface  charring  and  its  relation 
to  the  non-uniform  ablation  will  be  presented  in  a  parallel  paper 

The  plasma  density  and  electron  temperature  distribution  are  also  shown  in  Fig.  3.  The  plasma 
density  peaks  at  about  lO^"*  m"''  midway  between  the  electrodes.  The  electron  temperature  is 
strongly  non-uniform  radially  with  peak  near  the  electrodes  of  about  4.5  eV.  The  reason  for 
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higher  electron  temperature  near  the  electrodes  is  due  to  current  spreading  in  the  space  between 
the  electrodes  and  current  focusing  near  the  electrodes  (see  below  results  on  current  distribution). 
As  was  mentioned  earlier,  the  ablation  rate  is  also  non-imiform  radially.  This  effect  is  shown  in 
Fig.  4.  One  can  see  that  the  ablation  rate  peaks  near  the  electrodes  at  about  120  kg/m^s,  while  in 
the  middle  of  the  propellant  face  it  is  about  80-100  kg/m^s.  The  calculated  total  ablated  mass  per 
pulse  was  about  1  |Lig  that  is  close  to  the  measured  value  of  1.3  pg  [10]. 

A  region  of  magnetic  field  diffusion  in  the  near  field  outside  the  micro-PPT  is  shown  in  Fig.  5a. 
The  magnetic  field  drops  by  an  order  of  magnitude  at  about  1.5  mm  that  is  equal  to  the  thruster 
radius.  This  is  the  region  where  also  the  most  of  the  current  is  concentrated  as  shown  in  Fig.  5b. 
One  can  see  that  the  current  density  is  high  near  the  central  electrode  and  near  the  outer  electrode. 
This  is  a  reason  for  the  increasing  Teflon  surface  temperature  and  electron  temperature  in  these 
regions.  According  to  the  model  presented  in  Sec.  4  the  electromagnetic  acceleration  of  the 
plasma  is  also  expected  to  occur  in  this  region. 

Figure  6  shows  evolution  of  the  Carbon  ion  (C+)  component  of  the  plasma  plume  during  the  main 
part  of  the  pulse.  One  can  see  that  a  dense  plasma  focus  is  developed  at  a  few  mm  Ifom  the 
thmster  exit  plane.  This  plasma  focus  exists  during  the  entire  pulse  as  shown  in  Fig.  6,  but  the 
plasma  density  in  the  focus  decreases  from  about  2x10^^  m'^  at  the  beginning  of  the  pulse  down  to 
0.3x10^  m'^  at  5  ps.  At  the  beginning  (first  2  ps)  the  C+  density  mainly  develops  a  gradient  in 
the  radial  direction  that  is  a  result  of  high  directed  velocity  in  the  axial  direction.  Later,  during  the 
pulse,  the  axial  density  gradient  becomes  comparable  to  the  radial  one. 

The  Flourine  ions  (F+),  due  to  their  larger  mass,  have  different  dynamics  as  shown  in  Fig.7.  They 
have  smaller  acceleration  in  the  axial  direction  even  at  the  beginning  of  the  pulse  and  therefore 
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both  axial  and  radial  density  gradients  are  developed.  The  F+  density  in  the  plume  and  in  the 
plasma  focus  is  larger  than  that  of  C+,  because  originally  Teflon  has  composition  C2F4  with  F 
density  twice  larger  than  that  of  C.  Additionally  F  ions  experience  less  acceleration  in  the  plume 
because  of  their  mass  that  also  contribute  to  their  relative  density  increase. 

The  micro-PPT  is  essentially  an  electromagnetic  accelerator  as  shown  in  the  velocity  phase  plots 
(Figs.  8,9).  The  phase  plot  of  the  Carbon  ions  at  1  ps  is  centered  at  30  km/s  in  the  axial  direction. 
Ions  experience  also  radial  expansion  in  both  directions  due  to  the  magnetic  field  structure  and 
the  temperature  expansion.  The  radial  velocity  in  the  negative  direction  is  related  to  the  focus 
formation  along  the  axis,  as  shown  in  Figs.  6,7.  The  Flourine  ions  have  generally  smaller  both 
axial  and  radial  velocities  due  to  their  higher  mass.  At  a  later  stage  of  the  pulse  (see  Fig.  9) 
clearly  there  are  two  ion  populations  with  positive  and  negative  radial  velocities.  This  is  due  to 
the  annular  plasma  injection  corresponding  to  the  thruster  geometry  (see  Figs.  1,2). 

During  the  entire  pulse  there  is  a  population  of  ions  having  a  negative  axial  velocity  with  the 
magnitude  up  to  about  10  km/s  (see  Figs.  8,9).  This  population  creates  the  backflow 
contamination  that  is  an  important  issue  of  concern  for  a  spacecraft  using  the  PPT.  The  Carbon 
ions  have  a  larger  negative  velocity  due  to  their  higher  mobility  that  results  in  their  domination  in 
the  backflux. 
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7.  Summary 


In  this  paper,  a  self-consistent  description  of  an  electromagnetic  pulsed  plasma  thruster  from 
plasma  generation  into  the  near  plume  is  presented.  A  micro-PPT  developed  at  AFRL  is 
considered  as  a  working  example.  In  this  device,  no  separation  exists  between  the  main  plasma 
acceleration  region,  which  usually  occurs  in  an  internal  flow,  and  the  external  plasma  plume  field. 
Therefore,  a  single  end-to-end  model  is  necessary  for  accurate  simulations.  A  kinetic  Teflon 
ablation  model  is  incorporated  in  order  to  provide  the  boundary  conditions  for  the  plasma  plume. 
This  model  predicts  an  ablated  mass  per  pulse  of  about  1  p.g  that  is  close  to  that  measured  in 
experiment.  The  phenomena  in  the  plasma  plume  related  to  the  electromagnetic  effects  are 
studied.  The  plume  simulation  shows  that  a  dense  plasma  focus  is  developed  at  a  few  millimeters 
from  the  thruster  exit  plane  at  the  axis.  This  plasma  focus  exists  during  the  entire  pulse,  but  the 
plasma  density  in  the  focus  decreases  from  about  2x10^^  m'^  at  the  beginning  of  the  pulse  down 
0.3x10^^  m'^  at  5  p,s.  The  velocity  phase  is  centered  at  about  30  km/s  in  the  axial  direction 
demonstrating  that  the  micro-PPT  is  essentially  an  electromagnetic  accelerator.  At  a  later  stage  of 
the  pulse  there  are  two  ion  populations  with  positive  and  negative  radial  velocity.  It  is  predicted 
that  there  is  a  population  of  ions  having  a  negative  axial  velocity  magnitude  up  to  about  10  km/s. 
This  population  relates  to  the  backflow  contamination  that  is  an  important  issue  of  concern  for  a 
spacecraft  using  the  PPT. 
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Figure  3.  Teflon  surface  temperatxire,  plasma  density  and  elecrtron  temperature  distribution  in  the  layer 
near  the  Teflon  surface 
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Figure  6.  Evolution  of  the  Carbon  ion  density  during  the  pulse 
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Figure  9.  Ion  velocity  phase.  Late  stage  of  the  pulse 
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Abstract 

The  ablative  Z-pinch  PPT  utilizes  the  Z-pinch  effect  to  produce  an  axially  streaming  plasma.  When  the 
current  is  fully  pinched  in  this  device,  a  large  axial  pressure  gradient  exists  and  thus  plasma  accelerates  in 
the  axial  direction  due  to  the  gasdynamic  force.  In  the  present  paper,  a  model  of  the  electrical  discharge  in 
the  Ablative  Z-pinch  Pulsed  Plasma  Thruster  is  developed.  The  model  includes  Joule  heating  of  the  plasma, 
heat  transfer  to  the  Teflon,  and  Teflon  ablation.  Mechanisms  of  energy  transfer  from  the  plasma  colunrn  to 
the  propellant  include  heat  transfer  by  particle  convection  and  by  radiation.  The  computation  of  Teflon 
ablation  is  based  on  a  recently  developed  kinetic  ablation  model.  The  average  current  density  in  the  pinched 
region  is  used  as  a  parameter  of  the  model.  The  model  predicts  that  the  electron  temperature  peaks  at  about 
4  eV  and  the  plasma  density  peaks  at  about  8- 10^  m'^.  Thruster  performance  characteristics  such  as 
impulse  bit,  specific  impulse  and  the  mass  bit  are  calculated.  In  the  case  of  low  pulse  energy,  all  measured 
thruster  performance  characteristics  agree  with  our  model  predictions  when  the  average  current  density  to 
the  anode  current  density  ratio  (X  is  about  0.55.  In  the  case  of  high  pulse  energy,  such  an  agreement  with  the 
experiment  occurs  when  a=O.7-0.8,  that  suggests  that  in  the  case  of  higher  pulse  energy  the  current  is  more 
constricted  near  the  anode  tip.  The  model  also  predicts  that  the  impulse  bit  decreases  with  increasing 
propellant  inner  diameter  in  agreement  with  experiment.  The  comparison  of  the  model  prediction  with 
experimental  data  suggests  that  the  pinch  effect  and  the  thrust-to-power  ratio  increase  with  the  pulse 
energy. 
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1.  Introduction 


Due  to  the  combined  advantages  of  system  simplicity,  high  reliability,  low  average  electric  power 
requirement  and  high  specific  impulse,  currently  there  is  a  renewed  interest  in  pulsed  plasma 
thrusters  (PPT’s)  for  a  number  of  missions'.  The  PPT  is  considered  as  an  attractive  propulsion 
option  for  orbit  insertion,  constellation  maintenance,  drag  makeup  and  attitude  control  of  small 
satellites.  Existing  PPT’s,  however,  have  very  poor  performance  characteristics  with  an 

efficiency^  at  the  level  of  about  10%  that  leaves  an  opportunity  for  substantial  improvement. 

To  improve  the  PPT  performance,  several  directions  are  being  considered,  such  as  elimination  of 
late-time  ablation,  choice  of  the  proper  current  waveform  etcl  Currently,  new  PPT  devices  using 
both  an  electromagnetic  acceleration  mechanism''’^  or  an  electrothermal  mechanism  are  under 
development®’’.  One  of  the  motivations  for  development  of  new  PPT  configurations  is  to  achieve 
higher  thrust-to-power  ratio.  Electrothermal  PPTs  that  were  developed  by  Burton  at.  al.®’’  have 
thrust-to-power  ratio  >35  pNAV.  Another  approach  for  producing  a  high  thrust-to-power  PPT 
involves  using  the  inverse  Z-pinch  effect,  as  demonstrated  by  Mikellides  and  Turchi  et.  al.®. 

Z-pinch  plasmas  produced  by  the  magnetic  compression  of  a  cylindrical  plasma  column  are  used 
widely  to  produce  hot  and  dense  plasmas  for  many  applications  that  include  focusing  of 
energetic  particles,  guiding  intense  optical  laser  pulses,  and  producing  ultrahigh  magnetic  fields. 
While  the  macroscopic  dynamics  of  the  Z-pinch  plasma  can  be  described  in  terms  of  snowplow 
model  the  details  of  the  current  structure  and  other  plasma  properties  are  not  conq)letely 
understood’.  The  idea  to  use  the  Z-pinch  geometry  for  a  plasma  thruster  was  first  proposed  by 
Jahn  et.  al'°.  Recently  some  interesting  results  on  the  application  of  the  Z-pinch  configuration  for 
an  ablative  Pulsed  Plasma  Thruster  were  presented".  The  Ablative  Z-pinch  PPT  utilizes  the  pinch 
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effect  to  produce  axially  streaming  plasma.  It  was  found  that  a  specific  impulse  of  about  809  s,  a 
thmst-to-power  ratio  of  about  27.8  pNAV  and  a  thrust  efficiency  of  8.95%  were  the  highest 
performance  values  obtained  for  the  best  current  AZ-PPT  designs.  When  the  current  is  fully 
pinched  in  this  device,  a  large  axial  pressure  gradient  exists  and  thus  the  plasma  accelerates  in  the 
axial  direction.  Therefore,  the  main  mechanism  of  plasma  acceleration  is  electrothermal  due  to 
the  gasdynamic  force.  This  motivates  us  to  use  a  previously  developed  model  for  electrothermal 
PPT’s  to  describe  the  Z-pinch  PPT. 

Another  interesting  effect  that  may  enhance  the  AZ-PPT  performances  is  the  plasma 

macroparticle  interactions.  Macroparticles  are  large  chunks  (1-100  pm  diameter)  of  the  Teflon 

that  are  emitted  during  the  pulse.  Estimates  have  shown  that  the  particulate  emission  consumes 
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about  40%  of  the  total  propellant  mass,  while  contributing  only  1%  to  the  total  thrust  .  Specific 
design  of  the  AZ  PPT  allows  the  discharge  chamber  to  act  as  a  macroparticle  trap  as  was 
mentioned  in  Ref.  11.  Previously  we  showed  that  the  macroparticle  interaction  with  a  discharge 
plasma  may  lead  to  complete  decomposition  of  some  macroparticles  .  For  instance,  it  was 
found'^  that  a  5  pm  diameter  macroparticle  would  completely  decompose  during  the  10  ps  of  the 
interaction  with  a  plasma  under  typical  conditions  of  an  electrothermal  PPT.  Therefore 
macroparticle  trapping  can  increase  the  time  spent  by  the  macroparticle  in  the  discharge  that  will 
lead  to  increased  macroparticle  ablation  and  thus  better  propellant  utilization. 

In  this  paper,  we  describe  the  model  of  the  electrical  discharge  in  the  AZ-PPT.  Knowing  the 
plasma  parameter  evolution  during  the  pulse  allows  us  to  calculate  the  performance 
characteristics  of  the  thruster  such  as  specific  impulse,  impulse  bit  and  mass  bit.  The  present  work 
is  based  on  a  previously  developed  model  of  the  ablation  controlled  discharge''^'*^.  This  model 
was  successfully  used  to  model  elecrtrothermal  PPT’s  developed  at  the  University  of  Dlinois®’’' 


3 


(PPT-4,  PPT-7).  The  model  includes  Joule  heating  of  the  plasma,  heat  transfer  to  the  Teflon,  and 
Teflon  ablation.  Mechanisms  of  energy  transfer  from  the  plasma  column  to  the  propellant  bar 
include  heat  transfer  by  particle  convection  and  by  radiation.  The  computation  of  Teflon  ablation 
is  based  on  a  recently  developed  kinetic  ablation  model'®. 

In  the  following  section  we  will  describe  the  model  in  brief  that  includes  both  sheath  and  quasi¬ 
neutral  plasma  as  well  as  plasma-dielectric  interaction. 

2.  The  discharge  model 

The  model  considers  the  plasma  generation  processes  (ablation,  heating,  radiation,  ionization  etc.) 
and  plasma  acceleration  along  a  Teflon  chamber.  Some  characteristic  regions  such  as  the  Teflon 
surface,  electrical  sheath  near  the  dielectric  and  quasi-neutral  plasma  are  shown  in  Fig.  1. 
Different  kinetic  and  hydrodynamic  phenomena  determine  the  main  features  of  the  plasma  flow 
including  plasma  Joule  heating,  heat  transfer  to  the  dielectric  and  electrothermal  acceleration  of 
the  plasma  up  to  the  sound  speed  at  the  cavity  exit.  Below,  we  will  discuss  the  model  in  the 
different  regions  and  the  full  system  of  equations  including  the  final  expressions  obtained 
previously  [15].  In  addition,  we  will  estimate  the  current  distribution  in  the  pinched  area  in  order 
to  obtain  the  current  density  that  will  be  used  in  the  energy  balance. 

Firstly,  let  us  estimate  the  characteristic  time  for  the  pinch  effect.  In  the  frame  of  the  simple 
snowplow  model  when  the  current  rise  is  linear,  it  was  shown  that  the  time  constant  for  pinch'’: 

T  =  (fip)''^L/(dI/dt)’ 
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where  pi  is  the  permittivity,  p  is  the  initial  plasma  density,  dl/dt  is  the  current  rise,  and  L  is  the 
characteristic  initial  size.  For  the  range  of  parameters  typical  for  an  AZ-PFT,  an  estimation  shows 
that  the  current  pinch  time  is  <10  ^  s.  This  means  that  after  that  time,  the  current  is  strongly 
pinched.  This  statement  is  certainly  supported  by  the  experimental  observations  presented  earlier 
[  1 1].  Therefore,  in  our  model,  we  have  assumed  that  during  the  discharge  (~I0  ps)  the  current  is 
fully  pinched. 

2.1.  Electrostatic  Sheath 

According  to  our  previous  estimations*^  during  the  discharge  pulse,  a  quasi-steady  sheath 
structure  is  formed  and  that  under  typical  PPT  conditions,  this  sheath  is  unmagnetized  (the  self- 
magnetic  field  generated  during  the  pulse  is  considered).  In  this  case,  the  potential  drop  of  the 
electrostatic  sheath  near  the  Teflon  wall,  is  negative  in  order  to  repel  the  excess  of  the  thermal 
electrons,  so  that  the  random  electron  current  density  je*  is  equal  to  the  Bohm  ion  current  density 
ji.  Under  these  conditions  the  potential  drop  in  the  sheath  can  be  calculated  as; 


Ud  —  -Teln  (jetl/ji) 


(1) 


2.2.  Teflon  ablation 


In  the  present  work  the  Teflon  ablation  is  modeled  in  the  framework  of  the  approximation*®  based 
on  a  previously  developed  kinetic  model  of  metal  evaporation  in  a  surrounding  plasma  .  The 
mathematical  description  includes  the  model  for  two  different  layers  between  the  surface  and  the 
plasma  bulk:  (1)  a  kinetic  non-equilibrium  layer  adjusted  to  the  surface  with  a  thickness  of  about 
one  mean  free  path;  and  (2)  a  collision-dominated  layer  with  thermal  and  ionization  non- 


5 


equilibrium.  The  plasma-wall  transition  layer  also  includes  an  electrical  sheath  described  in  the 
previous  section.  This  model  makes  it  possible  to  calculate  the  plasma  parameters  (density  and 
temperature)  at  the  interface  between  the  kinetic  and  hydrodynamic  layers.  For  known  velocity 
and  density  at  this  interface,  it  is  possible  to  calculate  the  ablation  rate.  The  ablation  rate  is 
formulated  according  to  Ref.  16  as  follows: 

r  =  mV,n,  =  n,((2kT,/m)  (T2n2/2T,  -n,/2)/(nrn,%2)f ' .  (2) 

where  ni  and  Ti  are  the  density  and  temperature  at  the  kinetic  layer  edge,  and  n2,  T2  are  the 
density  and  temperature  at  the  hydrodynamic  layer-plasma  bulk  interface.  Density  na  and 
temperature  T2  are  determined  by  the  plasma  bulk  flow  and  energy  balance  (see  next  section). 
The  density  nj  and  temperature  Ti  are  determined  by  solution  of  the  problem  for  the  kinetic 
layer'*.  The  system  of  equations  is  closed  if  the  equilibrium  vapor  pressure  can  be  specified  that 
determines  the  parameters  at  the  Teflon  surface.  In  the  case  of  Teflon,  the  equilibrium  pressure 
formula  is  used  [1]: 

P  =  Pcexp(-Tc/Ts)=nskTs . (3) 

where  P  is  the  equilibrium  pressure.  Pc  =  1.84xl0“  N/m^  and  Tc  =  20815  K  are  the  characteristic 
pressure  and  temperature,  respectively. 

2.3.  Quasi-neutral  plasma  bulk 

The  schematic  geometry  of  the  thruster  is  shown  in  Fig.  1.  Several  simplifications  are  made  in 
order  to  make  it  possible  to  develop  a  simple  model.  For  instance,  the  anode  spike  is  assumed  to  a 
cylinder  as  shown  in  Fig.  1.  In  the  present  model  we  assume  that  all  parameters  vary  in  the  axial 
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direction,  x  (see  Fig.  1),  but  are  uniform  in  the  radial  direction.  The  axial  component  of  the  mass 
and  momentum  conservation  equations  reads: 


A(ap/at  +  a(pv)/3x) = 27iR2r(t,x) . w 

p(3V/0t  +V3V/ax)  =  -dP/dx . 


where  A  is  the  cross  section  of  the  Teflon  chamber  (A=7i(R2^-Ri^),  p  is  the  plasma  density,  P  is 
the  pressure,  V  is  the  plasma  velocity,  r(t,x)  is  the  local  instantaneous  ablation  rate,  Ra  is  the 
Teflon  chamber  radius  and  Ri  is  the  radius  of  the  spike  (see  Fig.  1). 

The  energy  transfer  from  the  plasma  column  to  the  wall  of  the  Teflon  cavity  consists  of  the  heat 
transfer  by  particle  fluxes  and  radiation  heat  transfer.  In  this  case,  the  energy  balance  equation 
can  be  written  in  the  form'^: 

|ne  OTe/at  +  vaTp/ax)  =  Qj  -  Qr  -  Qf . (6) 

where  Te  is  the  electron  temperature,  Qj  is  the  Joule  heat,  Qr  is  the  radiation  heat  and  Qf  is  the 
particle  flux.  The  radiation  energy  flux  Q,  includes  the  radiation  in  a  continuum  spectrum  based 
on  a  theoretical  model^°.  The  Joule  heat  source  is  assumed  to  be  concentrated  in  the  pinch  region 
by  the  abode  tip  (see  Fig.l).  The  average  current  density  in  this  region  is  used  as  a  parameter  of 
the  model.  The  particle  convection  flux  Qf  includes  energy  associated  with  electron  and  ion 
fluxes  to  the  dielectric  wall  that  leads  to  plasma  cooling.  Our  estimations  and  previous 
calculation  show^‘  that  the  electron  temperature  varies  only  slightly  with  axial  position  and 
therefore  we  performed  the  calculation  assuming  aTe/ax=0. 
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The  temperature  inside  the  Teflon  wall  can  be  calculated  from  the  heat  transfer  equation: 


8T/at  =  aa^/8r^, 


(7) 


where  a  is  the  thermal  diffusivity.  This  is  the  one-dimensional  equation  in  the  radial  direction. 
This  assumption  can  be  made  since  the  heat  layer  thickness  near  the  surface  is  smaller  than  the 
Teflon  cylinder  curvature  Rt  and  also  less  than  the  characteristic  length  of  plasma  parameter 
changes  in  the  axial  direction.  In  order  to  solve  this  equation,  we  use  the  following  boundary  and 
initial  conditions'^: 

-  A.aT/ax(x=0)  =  q(t)  -  AH-r  -  Cp(Ts-To)  r 

X8T/8x(x=oo)  =  0 . (g) 

T(t=0)=To 

where  x=0  corresponds  to  the  inner  dielectric  surface,  ZiH  is  the  ablation  heat,  F  is  the  rate  of 
Teflon  ablation  per  unit  area.  To  is  the  initial  room  temperature  and  q(t)  is  the  density  of  the  heat 
flux,  consisting  of  the  radiative  and  particle  convection  fluxes,  and  Tj  is  the  Teflon  surface 
temperature.  The  solution  of  this  equation  is  considered  for  two  limiting  cases  of  substantial  and 
small  ablation  rate  very  similar  to  that  described  in  Ref.  15. 

Having  calculated  the  plasma  density  and  electron  temperature,  we  calculate  the  chemical  plasma 
composition  considering  Local  Thermodynamic  Equilibrium  (LTE)  in  the  way  described 
previously'^'^^’^^.  In  the  considered  range  of  electron  temperature  (1-4  eV)  and  plasma  density 
(10“-10“  m'^)  we  will  assume  that  polyatomic  Teflon  molecules  C2F4  fully  dissociate  and  we 
will  start  our  consideration  from  the  point  when  we  have  gas  containing  C  and  F.  The  Saha 
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equations  for  each  species  (C  and  F)  are  supplemented  by  the  conservation  of  nuclei  and  quasi¬ 
neutrality. 

In  the  present  model  we  use  the  experimental  current  waveform  as  an  input  condition.  We  assume 
that  the  pulse  energy  is  supplied  by  a  simple  LRC  circuit  with  fixed  elements.  In  this  case,  the 
current  produced  by  an  underdamped  LRC  circuit  can  be  approximated  as: 


I(t)  =IpSin(at)exp(-Pt) 


(9) 


/2E  /  2 

where  Ip  =  A  / (x=  A  / p=  L  is  the  effective  inductance  in  the  circuit,  C  is  the 

capacitance,  R  is  the  total  circuit  resistance,  and  E  is  the  pulse  energy.  The  best  fit  with  the 
experimental  waveform  (frequency)  corresponds  to  {X=0.910*  s'*.  For  C=33.  6  pF  we  can 
estimate  that  L  in  the  circuit  is  about  3.7- 10'®  H.  For  an  energy  of  67  J,  the  amplitude  of  the 
current  waveform  Ip  equals  6- 10“*  A. 


The  current  waveform  was  measured  for  AZ  PPT-3  using  an  integrated  Rogowski  coil.  The  data 
that  were  taken  for  AZ-PPT3  fired  at  25, 50  and  67J  with  the  33.6  microF  capacitor  are  shown  in 
Fig.  2.  It  can  be  seen  that  the  thruster  exhibits  damped  sinusoidal  behavior  typical  of  ablative 
pulsed  plasma  thruster.  The  measurement  technique  was  the  same  as  described  in  Ref.  11.  The 
current  waveform  is  used  as  an  input  parameter  in  our  model. 


It  was  estimated  above  that  under  the  conditions  considered,  the  current  is  fully  pinched  during 
the  main  part  of  the  discharge.  The  current  is  concentrated  between  the  anode  tip  end  and  the 
cathode  as  schematically  shown  in  Fig.  1.  In  general,  the  current  density  distribution  in  that 
pinched  area  is  two-dimensional  and  can  be  calculated  using  the  magnetic  transport  equation. 
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However,  in  the  framework  of  the  present  ID  model  of  the  plasma  flow,  only  the  average  current 
density  is  considered.  The  current  density  peaks  at  the  anode  spike  tip,  ja=I/7tRi^  and  then 
decreases  toward  the  cathode.  In  the  present  work  we  use  the  current  density  as  a  parameter  in  the 
range  (0.2-  l)ja. 

3.  Results 

In  this  section  we  present  results  of  calculation  of  the  plasma  parameters  during  the  discharge 
pulse  for  the  AZ-PPT.  As  a  working  example,  we  consider  AZ-PPT-3  as  it  shows  the  best  device 
performance.  This  thruster  has  a  Teflon  chamber  length  of  57  mm,  the  radius  of  the  Teflon 
R2=12.5  or  19  mm,  and  the  radius  of  the  anode  spike  Ri~6  mm.  The  simulations  are  performed 
assuming  a  free  stream  condition  at  the  thruster  exit,  e.g.  the  plasma  velocity  equals  the  local 
sound  speed  at  x=L.  Based  on  the  calculated  plasma  parameter  distribution,  we  calculate  the 
thruster  performance  characteristics,  such  as  mass  ablated  during  the  pulse,  average  specific 
impulse,  and  gasdynamic  impulse  bit. 

The  plasma  density  spatial  and  temporal  distribution  is  shown  in  Fig.  3.  The  plasma  density  peaks 
at  about  8-10^^  m'^  at  3  ^is  that  corresponds  to  the  first  current  peak.  One  can  see  also  that  the 
plasma  density  decreases  towards  the  cathode  (along  x)  as  a  result  of  the  plasma  acceleration.  At 
the  exit  plane,  the  plasma  density  is  about  60%  of  the  plasma  density  near  the  anode  base  (x=0). 

The  electron  temperature  temporal  variation  is  shown  in  Fig.  4  with  discharge  energy  as  a 
parameter.  The  electron  temperature  peaks  at  3.7  eV  in  the  case  of  50  J  and  at  4.3  eV  in  the  case 
of  75  J  at  about  1.5  ps.  The  electron  temperature  oscillations  correspond  to  the  discharge  current 
oscillation  shown  above. 
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As  a  result  of  the  oscillations  in  electron  temperature  and  plasma  density,  the  ionization  degree 
calculated  assuming  LTE  also  oscillates  as  shown  in  Fig.  5.  Initially,  during  the  first  current  peak, 
the  plasma  is  strongly  ionized  and  the  thruster  produces  a  plume  containing  mainly  C  and  F  ions 
and  electrons.  At  late  time,  the  ionization  degree  at  the  peaks  is  about  0.5.  This  means  that  toward 
the  pulse  end  the  thruster  produces  a  large  amount  of  neutrals. 

The  calculated  thruster  performance  characteristics  integrated  over  the  15  fAs  pulse  are  shown  in 
Fig.  6.  Here  we  show  an  example  of  performance  calculation  for  AZ-PPT  3  firing  at  25  J.  For 
comparison  a  parameter  range  measured  experimentally  is  also  shown.  The  impulse  bit  and  mass 
bit  strongly  increase  with  parameter  a,  which  is  a  measure  of  the  current  density  in  the  pinched 
region.  One  can  see  that  all  measured  parameters  agree  with  our  model  predictions  when  the 
parameter  a=0.55.  This  result  suggests  that  even  though  the  current  is  focused  at  the  anode  tip, 
the  average  current  density  in  the  pinched  area  is  smaller  than  that  at  the  anode.  The  effect  of 
varying  the  Teflon  chamber  inner  diameter  (2R2,  see  Fig.  1)  is  shown  in  Fig.  7.  One  can  see  that 
the  model  predicts  a  decrease  of  impulse  bit  with  the  increasing  propellant  inner  diameter  (ID)  in 
agreement  with  experiment. 

The  calculation  of  the  AZ-PPT-3  performance  for  the  case  of  high  pulse  energy  is  shown  in  Fig.8. 
One  can  see  that  in  this  case,  agreement  with  the  experiment  occurs  when  (XM3.7-0.8.  This  means 
that  in  the  case  of  the  high  pulse  energy,  the  average  current  density  is  higher  and  closer  to  the 
anode  current  density.  This  is  an  expected  result  that  means  that  in  the  case  of  higher  pulse  energy 
(higher  current)  the  current  is  more  constricted  near  the  anode  tip.  The  effect  of  the  pulse  energy 
on  the  thrust-to-power  ratio  is  shown  in  Fig.  9  where  pulse  energy  is  used  as  a  parameter.  It  can 
be  seen  that  the  thrust-to-power  ratio  decreases  with  the  energy  for  the  constant  parameter  a. 
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However,  since  in  the  experiment  (Ref.  1 1)  it  was  obtained  that  the  thrust-to-power  ratio  actually 
increases  with  energy  one  can  conclude  that  the  parameter  a  must  be  higher  in  the  case  of  high 
energy.  Thus  the  comparison  of  the  model  prediction  with  experimental  data  suggests  that  the 
pinch  effect  increase  with  the  pulse  energy. 

The  simulated  results  are  summarized  in  Fig.  10  where  the  dependence  of  the  parameter  a 
(average  current  density  to  the  current  density  at  the  anode  tip  ratio)  for  which  the  agreement  with 
the  measured  impulse  bit  was  obtained  is  shown  as  a  function  of  the  pulse  energy.  Here  also  the 
Teflon  chamber  inner  diameter  (ID)  is  used  as  a  parameter.  Generally  it  can  be  seen  that  the 
average  cuirent  density  increases  with  the  pulse  energy  as  was  mentioned  above.  An  interesting 
result  is  the  dependence  of  the  average  current  density  on  the  Teflon  chamber  ID.  One  can  see 
that  the  average  current  density  more  significantly  increases  in  the  case  of  the  large  propellant  ID. 
This  observation  can  be  explained  as  follows.  When  the  propellant  ID  is  larger,  the  model 
predicts  that  the  plasma  density  in  the  channel  is  smaller.  Therefore  one  can  expect  higher  current 
density  in  the  case  of  larger  propellant  ID  since,  in  the  pinched  area,  the  current  constriction  is 
limited  by  the  plasma  pressure. 

Taking  into  account  that  the  average  current  density  varies  with  the  pulse  energy,  the  dependence 
of  the  thrust-to-power  ratio  (T/P)  on  the  Teflon  chamber  length  L  was  calculated.  These  results 
are  shown  in  Fig.  1 1 .  One  can  see  that  T/P  much  higher  in  the  case  of  the  smaller  propellant  ID. 
Also,  in  this  case  the  T/P  variation  with  the  energy  pulse  is  relatively  moderate.  In  the  case  of  the 
large  propellant  ID,  T/P  significantly  increases  with  the  pulse  energy.  In  all  cases  there  is  an 
optimal  length  corresponded  to  the  maximal  T/P.  In  the  case  of  1”  ID  propellant  the  maximal  T/P 
is  predicted  to  be  at  L-40  mm,  while  in  the  1.5”  ID  propellant  case  it  corresponds  to  L-50  mm. 
After  the  maximum  the  T/P  significantly  decreases  with  L,  especially  in  the  case  of  smaller  ED.  In 
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this  study  we  fix  the  length  of  the  pinched  area  meaning  that  when  the  Teflon  chamber  length  L 
increases  the  anode  spike  length  also  increases  (see  Fig.  1).  As  result  the  length  of  the  pinched 
area  normalized  by  the  total  length  of  the  propellant  decreases  and  therefore  the  plasma  in  the 
channel  is  less  heated.  This  is  the  main  reason  why  T/P  decreases  when  L  is  high. 

4.  Summary 

In  this  paper  a  model  of  the  discharge  in  the  recently  developed  Ablative  Z-pinch  Pulsed  Plasma 
Thruster  is  presented.  This  device  utilizes  the  pinch  effect  to  produce  an  axially  streaming 
plasma.  The  model  includes  Joule  heating  of  the  plasma,  heat  transfer  to  the  Teflon,  and  Teflon 
ablation  kinetics.  Mechanisms  of  energy  transfer  from  the  plasma  column  to  the  propellant  bar 
include  heat  transfer  by  particle  convection  and  by  radiation.  In  addition  it  was  assumed  that  the 
current  is  fully  pinched  near  the  anode  tip  during  the  main  part  of  the  discharge.  The  average 
current  density  in  the  pinched  area  is  used  as  a  parameter  of  the  problem.  The  model  predicts  that 
the  electron  temperature  peaks  at  about  4  eV  and  the  plasma  density  peaks  at  about  8-10  m  . 
During  the  initial  stage  of  the  discharge,  plasma  is  predicted  to  be  strongly  ionized,  while  at  the 
late  time  the  ionization  degree  peaks  at  about  0.5.  Thruster  performance  characteristics  such  as 
impulse  bit,  specific  impulse,  and  the  mass  bit  were  calculated.  In  the  case  of  low  pulse  energy, 
all  measured  thruster  performance  characteristics  agree  with  our  model  predictions  when  the 
average  current  density  to  anode  current  density  ratio  a  is  about  0.55.  In  the  case  of  high  pulse 
energy,  similar  agreement  with  experiment  occurs  when  a=0.7-0.8.  This  means  that  in  the  case  of 
higher  pulse  energy  (higher  current)  the  current  is  more  constricted  near  the  anode  tip.  The 
model  also  predicts  that  the  impulse  bit  decreases  with  increasing  propellant  inner  diameter  (ID) 
in  agreement  with  experiment.  The  comparison  of  the  model  prediction  with  experimental  data 
suggests  that  the  pinch  effect  and  the  thrust-to-power  ratio  increase  with  pulse  energy.  It  was 
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predicted  that  the  thrust-to-power  ratio  dependence  on  the  Teflon  chamber  length  has  a 
maximum.  The  non-monotonic  behavior  of  the  T/L  with  Teflon  chamber  length  is  more 
pronounced  in  the  case  of  the  smaller  propellant  inner  diameter. 


Acknowledgements 

The  authors  from  University  of  Michigan  gratefully  acknowledge  the  financial  support  by  the  Air 
Force  Office  of  Scientific  Research  through  grant  F49620-99- 1-0040. 


REFERENCES 

*  R.  L.  Burton  and  P.  J.  Turchi,  “Pulsed  plasma  thruster”,  Journal  of  Propulsion  and  Power, 
vol.14,  5, 1998,  pp.  716-735 

2r.  J.  Vondra,  ‘The  MIT  Lincoln  laboratory  pulsed  plasma  thruster”,  AIAA  Paper  76-998,  1976. 
^  P.  J.  Turchi,  Directions  for  improving  PPT  performance.  Proceeding  of  the  25'^  International 
Electric  Propulsion  Conference,  vol.  1,  Worthington,  OH,  1998,  pp.  251-258. 

*  E.  Antonsen,  R.L.  Burton  and  F.  Rysanek,  “Energy  measurements  in  a  co-axial  electromagnetic 
pulsed  plasma  thruster”.  Paper  AIAA- 1999-2292. 

^  F.  Gulczinski  III,  M.  Dulligan,  J.  Lakes  and  G.  Spanjers,  “Micropropulsion  research  at  AFRL”, 
Paper  AIAA-2000-3255 

®  R.L.  Burton  and  S.S.  Bushman,  “Probe  measurements  in  a  Co-axial  gasdynamic  PPT’,  35''' 
Joint  Propulsion  Conference,  Los  Angeles,  CA,  June  1999,  AIAA  Paper  99-2288. 

^  F.  Rysanek  and  R.L.  Burton,  “Effects  of  geometry  and  energy  on  a  coaxial  Teflon  pulsed 
plasma  thruster”,  36’’'  Joint  Propulsion  Conference,  Huntsville  AL,  July  2000,  AIAA  Paper  2000- 
3429. 

®  I.G.  Mikellides  and  P.J.  Turchi,  “Optimization  of  pulsed  plasma  thrusters  in  rectangular  and 
coaxial  geometries”,  EEPC  Paper  -  99-21 1. 

®  K.T.  Lee,  D.E.  Kim  and  S.H.  Kim,  “Reversed  current  structures  in  a  Z-pinch  plasma”.  Phys. 
Rev.  Lett.,  85,  2000,  3834. 

R.G.  Jahn,  W.  Jaskowsky,  and  R.L.  Burton,  “Ejection  of  a  pinched  plasma  from  an  axial 
orifice”,  AIAA  Journal,  3(10),  1965,  pp.  1862-1866. 


14 


"  T.E.  Markusic,  K.A.  Polzin,  J.Z.  Levine,  C.A.  McLeavey  and  E.Y.  Choueiri,  “Ablative  Z-Pinch 
Pulsed  Plasma  Thruster,  AIAA  Paper  2000-3257. 

12g.  G.  Spanjers,  J.  S.  Lotspeich,  K.A.  McFall,  and  R.  A.  Spores,  Propellant  losses  because  of 
particulate  emission  in  a  pulsed  plasma  thruster.  Journal  of  Propulsion  and  Power,  Vol.  14, 4, 
1998,  pp.  554-559. 

M.  Keidar,  I.  D.  Boyd  and  LI.  Beilis,  “Model  of  particulate  interaction  with  plasma  in  a  Teflon 
Pulsed  Plasma  Thruster”,  J.  Prop.  Power,  17,  No.  1,  2001,  pp.  125-131. 

M.  Keidar  and  I.D.  Boyd,  “Device  and  plume  model  of  an  electrothermal  pulsed  plasma 
thruster”.  Paper  AIAA-2000-3430. 

M.  Keidar,  I.D.  Boyd  and  LI.  Beilis,  “Electrical  discharge  in  the  Teflon  cavity  of  a  coaxial 
pulsed  plasma  thruster’,  IEEE  Trans.  Plasma  Sci.,  28,  2000,  p.  376-385. 

M.  Keidar,  LD.  Boyd  and  LI.  Beilis,  “On  the  model  of  Teflon  ablation  in  an  ablation-controlled 
discharge”,  J.  Phys.  D:  Appl.  Phys.,  34,  June,  2001,  pp.  1675-1677. 

"  N.  A.  Krall  and  A.W.  Trivelpiece,  Principles  of  Plasma  Physics,  McGraw-Hill,  New  York, 
1973. 

M.  Keidar,  J.  Fan,  LD.  Boyd  and  LI.  Beilis,  “Vaporization  of  heated  materials  into  discharge 
plasmas”,  J.  Appl.  Phys.,  89,  2001,  pp.  3095-3099. 

LI.  Beilis,  “Parameters  of  the  kinetic  layer  of  arc  discharge  cathode  region”,  IEEE  Trans. 
Plasma  Sci., 1985,  PS-13,  p.  288-290. 

“  G.  1.  Kozlov,  V.  A.  Kuznetsov,  and  V.  A.  Masyukov,  “Radiative  losses  by  argon  plasma  and 
the  emissive  model  of  a  continues  optical  discharge”,  Sov.  Phys.  JETP,  39, 1974,  pp.463-468. 

P.  Kovatya  and  J.  J.  Lowke,  ‘Theoretical  prediction  of  ablation  stabilised  arcs  confined  in 
cylindrical  tubes”,  J.  Phys.  D:  Appl.  Phys.,  17,  1984,  pp.  1197-1212. 

P.  Kovatya,  ‘Thermodynamic  and  transport  properties  of  ablated  vapors  of  PTFE,  alumina, 
perspex  and  PVC  in  the  temperature  range  5000-30000  K”,  IEEE  Trans.  Plasma  Sci.,  12,  1984 
pp.  38-42. 

C.S.  Schmahl  and  P.J.  Turchi, ,  “Development  of  equation-of-state  and  transport  properties  for 
molecular  plasmas  in  pulsed  plasma  thrusters.  Part  I:  A  two-temperature  equation  of  state  for 
Teflon”,  Proc.  Inter.  Electr.  Propul.  Conf.  Pp.  781-788, 1997. 


15 


Current  (A) 


Time, 


16 


0.1  0.2  0.3  0.4  0.5  0.6 

Axial  distance,  x/L 


Figure  3.  Plasma  density  temporal  and  axial  (along  channel  centerline)  variation  during  the  pulse. 
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Figure  6.  AZ-PPT-3  performances  as  a  junction  of  average  current  density  to  the  anode  current  density 
ratio  and  comparison  with  experiment 


Figure  7.  AZ-PPT-3  impulse  bit  as  a  junction  of  average  current  density  to  the  anode  current  density  ratio 
with  propellant  ID  (2R2)  as  a  parameter 
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Figure  8.  AZ-PPT-3  performances  as  a  function  of  average  current  density  to  the  anode  current  density 
ratio  and  comparison  with  experiment 


Figure  9.  ThrusMo-power  ratio  as  a  function  of  average  current  density  to  the  anode  current  density  ratio 
with  pulse  energy  as  a  parameter 
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Figure  10.  Average  current  density  to  the  anode  current  density  ratio  (X*  for  which  the  model  agree  with 
experiment  as  a  function  of  the  pulse  energy  with  capacitance  and  propellant  ID  as  parameters 
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Figure  II.  Thrust-to-power  ratio  as  a  function  of  the  Teflon  chamber  length  with  pulse  energy  and  the 
propellant  ID  as  parameters.  C=38.3  pF. 
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The  Teflon  ablation  in  a  micro-Pulsed  Plasma  Thruster  is  studied  with  an  aim  to 
understand  the  charring  phenomena.  Microscopic  analysis  of  the  charred  areas  shows  that 
it  contains  mainly  carbon.  It  is  concluded  that  the  carbon  char  is  formed  as  result  of 
carbon  flux  returned  from  the  plasma.  A  simplified  model  of  the  current  layer  near  the 
Teflon  surface  is  developed.  The  current  density  and  the  Teflon  surface  temperature  have 
peaks  near  the  electrodes  that  explain  preferential  ablation  of  these  areas  as  was  observed 
experimentally.  The  comparison  of  the  temperature  field  and  the  ablation  rate  distribution 
with  photographs  of  the  Teflon  surface  shows  that  the  area  with  minimum  surface 
temperature  and  ablation  rate  corresponds  to  the  charring  area.  This  suggests  that  the 
charring  may  be  related  to  a  temperature  effect.  Electron  densities  predicted  by  the  plume 
model  are  compared  with  near  field  measurements. 


Introduction 

Pulsed  plasma  thrusters  (PPT’s)  have  been 
investigated  since  the  early  1960’s  and  were  among 
the  first  of  various  electrical  propulsion  concepts 
accepted  for  space  flight  mainly  due  to  their  simplicity 
and  hence  high  reliability’.  However,  the  PPT  has  an 
efficiency  at  the  low  level  of  10%  (Ref.  2)  and 
therefore  several  ways  for  improvement  have  been 
suggested^.  Currently,  PPT’s  are  considered  as  an 
attractive  propulsion  option  for  stationkeeping  and 
drag  makeup  purposes  of  mass  and  power  limited 

4  5 

satellites  ’ .  Recently,  a  micro-PPT  has  being  designed 
at  the  Air  Force  Research  Laboratory  for  delivery  of 
very  small  impulse  bit^.  This  is  a  simplified 


miniaturized  version  of  a  conventional  PPT  with  a 
thrust  in  the  10  p,N  range  designed  to  provide  attitude 
control  and  stationkeeping  for  microsatellites. 

Complete  assessment  of  the  spacecraft  integration 
effects  requires  characterization  of  the  plasma  plume 
exhaust  of  a  PPT.  Previously  we  have  developed  an 
end-to-end  model  of  the  PPT  and  its  plume  with 
application  to  electrothermal^’^  and  electromagnetic 
PPT’s^.  It  became  clear  that  the  plasma  distribution  in 
the  plume  field  heavily  depends  upon  upstream 
boundary  conditions.  Therefore  the  model  of  the 
plasma  generation  in  these  devices  becomes  a  very 
important  aspect  of  accurate  plasma  plume  simulation. 


Copyright  ©  2001  by  Michael  Keidar.  Published  by  Electric  Rocket  Propulsion  Society  with  permission 


1 


In  the  present  paper  we  will  focus  on  the  MicroPPT. 
Inspection  of  the  micro-PPT  propellant  surface  after 
firing  indicated  signs  of  charring  and  preferential 
ablation  near  the  electrodes*’'®.  We  present  also  results 
of  the  microscopic  analyses  of  the  charring  areas  that 
is  the  useful  tool  for  understanding  of  the  charring 
mechanism.  In  order  to  understand  this  phenomenon  a 
model  of  the  plasma  layer  near  the  Teflon  surface  is 
developed.  In  addition,  the  solution  of  the  model  will 
provide  boundary  conditions  for  the  plasma  plume. 

Microscopic  analyses  of  the  Teflon  surface 

The  AFRL  MicroPPT  in  development  for  TechSat21 
uses  a  3  electrode  configuration*.  A  small  diameter 
rod  (center  electrode)  is  encased  in  a  small-diameter 
annulus  of  Teflon,  which  is  then  encased  in  a 
relatively  small  diameter  tube,  which  acts  as  the 
intermediate  electrode.  This  construction  is  then 
encased  in  a  second  larger  diameter  annulus  of  Teflon, 
which  is  then  encased  in  a  large  diameter  outer 
electrode.  The  MicroPPT  is  fired  by  a  low-energy 
breakdown  between  the  intermediate  and  central 
electrode.  This  discharge  provides  enough  seed 
ionization  to  enable  the  higher  energy  conduction 
breakdown  between  the  intermediate  and  out 
electrodes.  The  discharge  between  the  intermediate 
and  central  electrode  is  referred  to  as  the  trigger 
discharge.  The  discharge  between  the  intermediate 
and  outer  electrode  is  referred  to  as  the  “main 
discharge”.  Although  a  wide  range  of  parameters  are 
tested  in  various  MicroPPT  configurations,  typically 
the  trigger  discharge  will  consiune  about  1/50  the 
energy  of  main  discharge.  In  this  fashion,  the 
MicroPPT  has  demonstrated  the  ability  to  passively 
initiated  a  surface  breakdown  discharge  across  outer 
propellant  diameters  as  high  as  14”  using  a  relatively 
low  voltage  below  3000V.  Without  the  3-electrode 
configuration,  up  to  40  kV  would  be  required  to 
initiate  the  discharge  across  a  14”  diameter.  Requiring 
a  40  kV  charge  would  place  excessive  design 
requirements  on  the  power-processing  unit  and  on  the 
spacecraft  EMI  shielding. 

In  this  work,  research  is  performed  on  2-electrode 
MicroPPT  configurations.  The  discharge  occurs 
between  an  inner  cathode  rod  and  an  outer  anode  tube, 
across  a  Teflon  annulus.  Understanding  the  physical 
processes  in  this  simplified  geometry  has  proven 
beneficial  in  advancing  the  optimization  of  the 
MicroPPT  by  separating  the  requirements  for  the 


trigger  and  main  discharges.  Research*’’®  on  small 
diameter  2-electrode  designs,  generally  between  1-3 
mm,  is  applicable  to  the  trigger  discharge. 
Research  *’"  on  larger  diameter  2-electrode  designs, 
typically  between  3  and  7  mm,  are  more  applicable  to 
the  main  discharge. 

Micro-PPT  propellant  samples  with  2  electrodes  and 
different  anode  diameters  were  analyzed.  These 
samples  represent  a  fully  charred,  a  partially  charred 
and  an  imcharred  Teflon  surface.  Microscopic 
analyses  were  performed  on  the  Environmental 
Scanning  Electron  Microscope  available  at  the  EMAL 
Center  at  the  University  of  Michigan.  The  sample 
chamber  is  held  at  a  pressure  of  typically  between  1- 
20  torr.  The  accelerating  voltage  is  about  15-20  kV. 
X-ray  Energy  Dispersive  Spectroscopy  (XEDS) 
makes  it  possible  to  identify  the  chemical  elements. 
Below  we  present  some  characteristic  images  taken 
from  the  fully  charred  sample. 


Figure  L  Charring  area  on  the  propellant  surface  of  3.6 
mm  micro-PPT 


Generally  two  different  structures  are  identified  in  the 
charred  area  (shown  as  Area  1  and  Area  2  in  Fig.  1)  as 
can  be  concluded  also  from  the  high  magnification 
images.  One  can  conclude  that  Area  1  mainly  contains 
Carbon  and  other  small  fractions  of  Fluorine,  Copper 
(Cu),  Silver  (Ag)  and  Silicon  (Si)  as  shown  in  Fig.  2. 
The  images  from  Area  2  look  very  different  as  shown 
in  Figure.  3.  One  can  see  that  the  main  component 
here  is  Silicon.  Small  fractions  of  Cu  and  Ag  are  also 
found. 
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An  analysis  of  the  interface  between  Areas  1  and  2  as 
shown  in  Figure  1  suggests  that  the  same  structure  (as 
in  the  image,  Fig.  3)  may  lie  under  the  Carbon 
charring.  In  order  to  verify  this  we  removed  the 
Carbon  layer  and  analyzed  the  scratched  area.  We 
conclude  from  the  element  mapping  that  the  scratched 
area  contains  Silicon  as  well  as  Fluorine. 


Figure  2.  Image  and  XEDS  results  for  the  Area  L  The  main  peak 
corresponds  to  Carbon 

Microscopic  analyses  of  different  fully  and  partially 
charred  samples  show  that  under  the  Carbon  layer 
there  is  a  layer  of  Silicon  with  some  small  amount  of 
Copper.  The  origin  of  Copper  is  probably  the  outer 
electrode  while  Silicon  may  come  from  the  diffusion 
pump  or  the  vacuum  facility. 

In  order  to  eliminate  this  possible  source  of  Silicon  the 
MicroPPT  was  fired  in  a  chamber  with  a  turbopump 
(glass  bell  jar).  The  image  of  this  sample  is  shown  in 
Fig.  4  and  the  typical  image  of  the  charred  area  is 
shown  in  Fig.  5. 


Figure  3.  Image  and  XEDS  results  for  the  Area  2.  The  main  peak 
corresponds  to  Silicon 


Figure  4.  Propellant  surface  of  3. 6  mm  diameter  micro- 
PPT 

It  is  interesting  to  note  that  in  the  micro-PPT  sample 
fired  in  the  turbopumped  chamber,  there  is  no 
evidence  of  Silicon. 
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An  important  observation  from  the  microscopic 
analyses  is  the  presence  in  most  samples  of  a  layer  of 
metal  under  the  char.  In  those  cases  where  no  metal 
layer  is  found  under  the  char,  the  charred  area  has  the 
same  appearance.  It  is  concluded  that  the  char 
formation  therefore  may  be  the  same  in  both  cases. 
This  fact  may  suggest  that  the  Carbon  char  is  formed 
as  result  of  the  Carbon  flux  returned  from  the  plasma 
rather  than  non-complete  Teflon  decomposition. 


Figure  5.  Charred  area  image  and XEDS  results.  The  main  peaks 
corresponds  to  Carbon  and  Copper 


Model 

In  this  part  we  describe  a  model  of  the  plasma  layer 
near  the  evaporating  surface  with  application  to  a 
micro-PPT  that  is  shown  schematically  in  Fig.  6. 

The  model  includes  the  following  features;  Teflon 
ablation,  plasma  energy  balance,  heat  transfer  from 
the  plasma  to  the  Teflon,  current  spreading  in  the  near 
field,  and  an  equivalent  RLC  electrical  circuit  model. 
The  Teflon  ablation  model  is  based  on  a  recently 
developed  kinetic  ablation  model.'^  '^ 


Figure  6.  Schematic  of  the  problem  geometry  and  energy  balance 

The  plasma  energy  balance  and  the  heat  transfer  to  the 
Teflon  are  based  on  a  previously  developed  model  of 
the  ablation  controlled  discharge.*’’'*  The  energy 
balance  in  the  layer  has  the  following  form: 

3/2ndT/dt=  Qj-Qr-Qr . (1) 

where  n  is  the  plasma  density,  T  is  the  plasma 
temperature,  Qj  is  the  Joule  heat,  Qr  is  the  radiation 
heat  and  Qf  is  the  energy  flux  due  to  particle 
convection.  Mechanisms  of  energy  transfer  from  the 
plasma  column  to  the  propellant  bar  include  heat 
transfer  by  particle  convection  and  by  radiation.  The 
inputs  for  the  model  are  thruster  geometry.  Teflon 
material  properties,  and  Teflon  equilibrium  pressure 
dependence  on  the  surface  temperature. 

In  the  transition  region  between  the  plasma  and  the 
ablated  surface,  two  different  layers  are  distinguished: 
a  kinetic  non-equilibrium  layer  adjusted  to  the  surface 
with  a  thickness  of  about  one  mean  free  path;  and  (2) 
a  collision-dominated  layer  with  thermal  and 
ionization  non-equilibrium.  The  solution  for  these  two 
layers  is  coupled  with  the  quasi-neutral  plasma  that 
allows  the  calculation  of  the  ablation  rate.  The  energy 
input  in  Eq.  1  depends  upon  the  current  density 
distribution  near  the  ablated  surface.  In  order  to 
calculate  the  current  density,  the  problem  of  the 
current  distribution  in  the  thruster  near  field  is  solved. 

Current  distribution  in  the  plasma  plume  near 
field 

Assuming  that  the  magnetic  field  has  only  an 
azimuthal  component  and  after  neglecting  the 
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displacement  current,  the  magnetic  field  in  the  near 
field  plasma  plume  is  calculated  from  the  magnetic 
transport  equation  in  the  following  form: 

aB/at  =  l/(a|i)V^B  -  VxO'xB/en)  +  Vx(VxB)  (2) 

where  a  is  the  plasma  conductivity,  |ii  is  the 
permittivity,  j  is  the  current  density  and  V  is  the 
plasma  velocity. 

A  scaling  analysis  shows  that  the  various  terms  on  the 
right  hand  side  of  Eq.  2  may  have  importance  in 
different  regions  of  the  plasma  plume  and  therefore  a 
general  end-to-end  plasma  plume  analysis  requires 
keeping  all  terms  in  the  equation.  In  the  case  of  the 
near  plume  of  the  MicroPPT  with  a  characteristic 
scale  length  of  about  1  cm,  the  magnetic  Reynolds 
number  Rem«l  and  therefore  the  last  term  can  be 
neglected.  Taking  this  into  account  in  the 
dimensionless  form,  Eq.  2  can  be  written  as: 

RemdB/at  =  V^B  -  (tOT)  {Vx(VxBxB)}  (3) 


here  (cor)  is  the  Hall  parameter  that  measures  the  Hall 
effect.  Therefore,  depending  on  the  plasma  density, 
the  Hall  effect  may  be  important  for  the  magnetic  field 
evolution. 


Figure  7.  Magnetic  field  distribution  in  the  near  field  of  a 
3. 6  mm  diameter  micro-PPT 
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Our  estimations  show  that  the  Hall  parameter  cox«l 
if  the  plasma  density  near  the  Teflon  surface 
N>10^^  m‘^.  However,  in  general,  the  Hall  parameter 
may  vary  for  different  devices  and  therefore  in  the 
future  we  will  investigate  the  effect  of  the  Hall 
parameter  on  the  magnetic  field  distribution  in  the 
near  field.  From  the  magnetic  field  distribution  the 
current  density  components  can  be  calculated  as 
follows: 

Jr  =  -l/p9B/9z 
Jz=l/|Li  dB/dr 

An  example  of  the  magnetic  field  distribution  is 
shown  in  Fig.  7  for  a  3.6  mm  diameter  propellant 
micro-PPT.  The  current  density  distribution  is  shown 
in  Fig.  8.  One  can  see  that  the  thickness  of  the  layer 
where  the  main  part  of  the  current  is  concentrated  is 
about  1  mm  from  the  Teflon  surface.  More  details 
about  the  current  and  magnetic  field  distributions  in 
the  near  field  are  presented  in  a  recent  paper^. 


Figure  8.  Current  distribution  in  the  near  field  of  a  3.6  mm 
diameter  micro-PPT 

The  current  density  radial  distribution  (according  to 
Eqs.  2-4)  near  the  Teflon  surface  is  shown  in  Fig.  9. 


Radial  distance,  mm 

Figure  9.  Current  density  near  the  Tefion  surface 

One  can  see  that  in  the  case  of  the  smaller  thruster 
(2.21  mm  diameter)  the  current  density  is  higher  by  a 
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factor  of  about  2  for  the  same  total  discharge  current. 
The  current  density  near  the  Teflon  surface  has  a 
minimum  due  to  the  current  spreading  in  r-z  plane  as 
shown  in  Fig.  8. 

The  dependence  of  the  calculated  ablated  mass  during 
the  pulse  is  shown  in  Fig.  10  for  two  micro-PPT’s 
with  diameters  of  3.6  mm  and  2.21  mm. 


Figure  10.  Ablated  mass  as  a  function  of  time  with 
propellant  diameter  as  a  parameter 

It  should  be  noted  that  in  the  experiments,  the  average 
ablation  mass  was  measured  to  be  about  1.3  ^g/shot  in 
the  case  of  the  2.21  mm  diameter  thruster.  Our  model 
predicts  without  any  fitting  parameter  that  the  mass 
ablated  per  pulse  is  about  0.9  p,g  in  this  case.  Taking 
into  account  that  the  late  ablation  in  the  form  of 
macroparticles  (that  is  not  considered  here  but  was 
observed  in  the  micro-PPT*^  may  consume  up  to  40% 
of  the  mass*^,  one  can  conclude  that  the  model 
predicts  the  ablation  rate  reasonably  well. 

According  to  the  energy  balance  (Eq.  1)  and  the  heat 
transfer  equation  at  the  Teflon  surface,  the  surface 
temperature  should  depend  on  the  current  density.  The 
spatial  and  temporal  variation  of  the  Teflon  surface 
temperature  for  the  two  thrusters  is  shown  in  Fig.  1 1 . 
In  these  calculations,  the  experimental  current 
waveform  is  used.  In  the  case  of  the  thruster  with 
smaller  diameter  (2.21  mm),  the  Teflon  surface 
temperature  is  higher  by  about  20  K  and  can  be 
considered  more  uniform  radially  than  that  of  the 
thruster  with  larger  diameter  during  the  whole  pulse. 
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Figure  11.  Teflon  surface  temperature  (K)  temporal  and  radial 
variation  for  two  MicroPPT  designs  with  diameter  2,21  mm  and 
3.6  mm 

In  the  larger  thruster  (3.6  mm),  the  temperature  has  a 
minimum  at  radial  distances  of  1.1 -1.3  mm.  Since  the 
Teflon  ablation  is  approximately  exponentially 
proportional  to  the  surface  temperature,  the  model 
predicts  a  lower  rate  of  ablation  in  the  areas  where  the 
surface  temperature  has  a  minimum.  Taking  this  into 
account,  the  effect  of  the  temperature  distribution  may 
be  related  to  the  preferential  charring  of  the  Teflon 
surface  observed  experimentally  as  shown  in  Fig.  12. 
It  is  interesting  to  note  that  comparison  of  the 
calculated  temperature  field  and  ablation  rate  with  the 
photograph  of  the  Teflon  surface  (see  Fig.  12)  shows 
that  the  area  with  surface  temperature  and  ablation  rate 
minimum  corresponds  to  the  charring  area  in  the  case 
of  the  3.6  mm  diameter  thruster. 
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Figure  12,  Teflon  surface  photo  (Ref.  6)  and  the  Teflon  surface 
temperature  field  and  the  ablation  rate  in  the  case  of  3.6  mm 
diameter  micro-PPT 


Effect  of  discharge  energy 


The  current  produced  by  an  underdamped  RLC  circuit 
has  the  following  form: 

I(t)  =  Ip  sin(at)-exp(-pt) 

where  a=(l/LC  -  P^),  P=R/2L  and  Ip  is  the  current 
peak,  Ip=Uo/(La),  where  Uo  is  the  initial  voltage  on 
the  capacitor  of  capacitance  C,  R  is  the  equivalent 
circuit  resistance,  L  is  the  circuit  inductance.  From  the 
comparison  of  I(t)  with  the  experimental  current 
waveform  in  the  case  of  C=0.3  |uiF  it  is  estimated  that 
R=0.3  Q  and  L=3.610'’  H. 

We  studied  the  effect  of  the  discharge  energy  (CUo^) 
on  the  Teflon  surface  temperature  and  the  Teflon 
ablation  rate.  These  results  are  shown  in  Fig.  13.  One 
can  see  that  with  the  energy  increase,  the  Teflon 
surface  temperature  and  the  ablation  rate  increase. 
These  results  suggest  that  increase  of  the  discharge 
energy  for  constant  capacitance  leads  to  enhanced 
Teflon  ablation. 
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Figure  13,  Ablation  rate  temporal  and  radial  variation  for 
discharge  energies  of  3J  and  6J. 

It  was  shown  in  experiment  that  the  energy  level 
affects  the  Teflon  ablation.  These  experiments  were 
conducted  in  the  glass  bell  jar  under  a  pressure  of  10'^ 
torr.  The  char  patterns  for  three  different  energies  are 
shown  in  Fig.  14  after  continuous  firing  for  at  least  8 
hours  and  at  1  Hz.  The  cases  with  higher  discharge 
energy  were  fired  for  longer  duration  to  see  if  char 
would  appear.  There  have  been  no  experimental 
observations  of  cases  where  char  appeared,  only  to  be 
cleaned  up  by  subsequent  discharges  at  the  same 
energy.  (Cleaning  of  the  char  formation  by  firing  at 
higher  discharge  energies  has  been  observed 
experimentally,  but  is  not  considered  within  the 
context  of  this  effort).  For  these  tests  a  diameter 
2-electrode  MicroPPT  was  energized  using  a  0.417  uF 
capacitor.  The  charge  voltage  ranged  from  2448V  for 
the  1.25  J  case,  to  5364V  for  the  6  J  case.  Clearly 
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these  voltages  are  insufficient  to  cause  the  surface 
breakdown  needed  for  MicroPPT  discharge  initiation 
across  a  14”  diameter.  Rather  than  complicate  the  test 
setup  by  using  the  3 -electrode  MicroPPT 
configuration*,  an  auxiliary  sparkplug,  fired  at  0.5  J, 
was  used  to  initiate  the  discharge  on  command. 

One  can  see  that  in  the  case  of  small  energy,  the 
charring  in  the  area  between  the  electrodes  is  observed 
while  in  the  case  of  large  energy,  there  is  no  charring. 
The  middle  energy  case  shows  a  level  of  char  between 
the  2  extrema.  This  effect  can  be  partially  explained 
in  terms  of  om  model  that  shows  that  higher  discharge 
energy  leads  to  higher  Teflon  surface  temperature 


Figure  14,  Photo  of  the  micro-PPT  propellant  surface  for 
discharge  energies  of  1.25  J,  3J  and  6J. 


Near  field  plume 

In  this  section  we  will  present  measured  and  predicted 
electron  density  distributions  in  the  near  filed  plume 
for  one  micro-PPT  design.  These  data  will  be 
compared  in  order  to  verify  our  plume  and  device 
model. 

Herriot  Cell  electron  density  measurement 

An  experimental  basis  for  comparison  is  provided 
using  a  Herriott  Cell  interferometer.  Electron  density 
measurements  are  taken  on  a  6.35  mm  (1/4”)  diameter 
MicroPPT  at  AFRL.  The  interferometer  uses  a  single 
laser  wavelength  and  quadrature  heterodyne  technique 
described  by  Spanjers  et  al 

Addition  of  a  Herriott  Cell  acts  to  confine  a  large 
number  of  laser  passes  into  an  area  suitable  for 
maximum  exposure  to  the  MicroPPT  plume.  This  is 
achieved  by  focusing  the  laser  between  the  two 
concave  mirrors  of  the  cell.  The  technique  is  used  to 


increase  signal-to-noise  ratio  for  diffuse  plasmas  by 
increasing  laser  exposure  to  the  plasma  over  a 
characteristic  path  length.*^  Thirteen  laser  reflections 
in  the  Herriott  Cell  were  focused  to  two  points, 
separated  by  3  mm.  For  data  shown  here,  these  points 
formed  a  plane  parallel  to  the  fuel  face  and  5  mm 
distant.  A  schematic  of  the  beam  geometry  is  shown 
in  Fig.  15. 


Figure  15.  Schematic  showing  interferometer  coverage  of  the 
Micro-PPT  fuel  face  using  the  Herriott  Cell.  The  beams  are 
situated  5mm  from  the  propellant  face. 

Figure  16  shows  the  experimental  data  from  this 
geometry  co-plotted  with  model  predictions.  The 
experimental  data  was  taken  at  a  discharge  energy  of 
6,6  J  from  a  0.417  pF  capacitor.  Experimental 
waveforms  of  the  current  were  obtained  using  a  a  self- 
integrating  Rogowski  coil.  Peak  density  reaches 
23±6xl0*^  cm"^  with  uncertainty  due  to  shot-to-shot 
variations  in  thruster  firing. 


Near  field  plume  simulations 


Using  the  plasma  layer  model  predictions  as  boundary 
conditions,  we  calculated  the  near  field  plume  of  the 
MicroPPT.  This  allows  us  to  make  direct  comparison 
of  our  model  predictions  with  measured  data. 

The  general  approach  for  the  plume  model  is  based  on 
a  hybrid  fluid-particle  approach  that  was  used 
previously  (Refs.  7).  In  this  model,  the  neutrals  and 
ions  are  modeled  as  particles  while  electrons  are 
treated  as  a  fluid.  Elastic  (momentum  transfer)  and 
non-elastic  (charge  exchange)  collisions  are  included 
in  the  model.  The  particle  collisions  are  calculated 
using  the  direct  simulation  Monte  Carlo  (DSMC) 
method^^.  Acceleration  of  the  charged  particles  is 
computed  using  the  Particle-In-Cell  method  (PIC)^^. 
The  ion  dynamics  is  calculated  by  taking  into  account 
electromagnetic  acceleration. 
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The  electron  dynamics  is  very  important  in  the  plasma 
plume  of  an  electromagnetic  PPT.  Previously,  our 
model  was  based  on  the  assumption  that  electrons 
rapidly  reach  the  equilibrium  distribution  and  in  the 
absence  of  a  magnetic  field  can  be  described 
according  to  the  Boltzmann  distribution.  While  this 
was  a  satisfactory  assumption  in  the  case  of  an 
electrothermal  thruster  plume, ^  this  is  not  suitable  for 
the  near  field  of  an  electromagnetic  thruster.  In  the 
case  of  a  magnetic  field,  the  electron  momentum 
equation  reads  (neglecting  electron  inertia): 

0  =  -e^ne(E+VeXB)  -  eVPe  -  Veimj 

The  electric  and  magnetic  field  distributions  in  the 
plume  are  calculated  from  the  set  of  Maxwell 
equations.  More  detailed  study  of  the  near  field  plume 
of  the  electromagnetic  PPT  was  presented  recently^. 


Figure  16.  Comparison  of  predicted  and  measured  electron 
density  time  variation  at  5  mm  from  the  propellant  face  at  the  axis 
in  the  case  of  6.35  mm  diameter  micro-PPT firing  at  6.6  J. 

The  calculated  plasma  density  time  variation  at  5  mm 
from  the  propellant  face  at  the  axis  is  shown  in  Fig. 
16.  Plasma  density  peaks  at  about  3x10^^  m'^  and 
decreases  by  few  order  of  magnitude  towards  the 
pulse  end.  For  comparison  we  show  also  measured 
plasma  density  (using  Herriott  Cell  technique,  see 
above  section).  One  can  see  that  our  model  correctly 


predicts  both  the  plasma  density  level  and  temporal 
behavior  during  the  entire  pulse. 

Summary 

A  model  of  the  plasma  layer  near  the  Teflon  surface 
of  a  PPT  was  developed  that  allows  the  calculation  of 
the  Teflon  surface  temperature  and  ablation  rate  self- 
consistently.  It  was  found  that  the  propellant  size  has 
an  important  effect  on  the  Teflon  surface  temperature 
distribution  and  the  ablation  rate.  For  instance,  in  the 
case  of  the  thruster  with  smaller  diameter  (2.21  mm) 
the  Teflon  surface  temperature  is  higher  by  about  20 
K  and  can  be  considered  more  uniform  radially  than 
that  of  the  thruster  with  larger  diameter  during  the 
whole  pulse,  hi  the  larger  thruster  (3.6  mm),  the 
temperature  has  a  minimum  at  radial  distances  of  1.1- 
1.3  mm.  The  comparison  of  the  temperature  field  and 
the  ablation  rate  distribution  with  a  photograph  of  the 
Teflon  surface  shows  that  the  area  with  surface 
temperature  and  ablation  rate  minimum  corresponds 
to  the  charring  area  in  the  case  of  the  3.6  mm  thruster. 
This  suggests  that  the  charring  may  be  related  to  the 
temperature  effect.  An  analysis  of  the  effect  of  the 
discharge  energy  E  on  the  temperature  distribution 
shows  that  the  Teflon  surface  temperature  and  the 
ablation  rate  can  be  increased  by  increasing  E.  At  the 
same  time,  the  increase  of  capacitance  leads  generally 
to  a  smaller  ablation  rate,  though  this  effect  can  be 
considered  to  be  marginal. 

A  microscopic  analysis  of  the  charred  areas  showed 
that  the  charred  area  contained  mainly  Carbon.  In 
some  cases  a  metal  layer  was  found  under  the  Carbon 
char.  The  metal  deposition  is  related  to  the  electrode 
erosion  while  the  Silicon  is  assumed  to  come  from  the 
diffusion  pump.  In  fact,  when  a  cryogen  pumping 
system  was  used,  no  Silicon  was  obtained  on  the 
Teflon  surface.  In  those  cases  where  no  metal  layer 
was  found  under  the  char,  the  charred  area  has  the 
same  appearance.  It  is  concluded  that  the  char 
formation  therefore  may  be  the  same  in  both  cases. 
This  fact  may  suggest  that  the  Carbon  char  is  formed 
as  result  of  the  Carbon  flux  returned  from  the  plasma 
rather  than  non-complete  decomposition  of  the  Teflon. 

Predicted  electron  density  was  directly  compared  with 
experimental  data  and  very  good  agreement  was 
obtained. 
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